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CHAPTER I 


INTRODUCTION 

The purpose of this thesis is the presentation of an improved 
method for obtaining numerical solutions of a certain class of two-point 
boundary value problems which often arise in optimal control theory. 

These problems are characterized by systems of nonlinear ordinary differ- 
ential equations with nonlinear boundary conditions. 

A general problem in optimal control theory is often stated in the 
following manner. Given a system which is described by a set of non- 
linear ordinary differential equations 

. x = f(x,u,t) (l.l) 

where x is an n vector describing the state (position and velocity) 
of the system as a function of time t , u is a q vector of time 
varying controls which can be applied to the system, and f is an 
n vector of nonlinear functions; it is required to determine the control 
histories u(t), so that an extremum (maximum or minimum) of some scalar 
performance index <f>( x(t^) ,t f ) is obtained at some terminal time t f . 

The system must satisfy certain initial conditions in the form of an 
m vector of functions 

-5-( x -(-o) ’- o) -- ° - - - - >- 
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and the controls must not only extremize the performance index, but also 
yield a final system state which satisfies a k vector of nonlinear 
functions 

ip (x ,t f ^ =0 k < n (1.3) 

The problem as stated is known as a Mayer problem in the calculus of 
variations. Simple transformations are given by Bliss [l] which trans- 
form the Lagrange problem (where the performance index is a definite 
integral) and the more general Bolza problem into Mayer problems. With 
these transformations, a large number of optimal control problems can 
be stated in the convenient Mayer form. 

Since the control problem stated is rarely amenable to analytic 
solution, methods for obtaining numerical solutions are necessary. The 
availability of large digital computers coupled with a demand for solu- 
tions in various space age applications has resulted in considerable 
research activity in the solution of the optimal control problem by 
numerical methods. The majority of this work has occurred in the past 
ten years and methods of solution are basically divided into either 
direct or indirect methods. 

The direct methods are so-called because they seek to directly 
manipulate the control histories u(t) in order that an augmented 
functional (which includes the given performance index and a measure 
of terminal constraint satisfaction) is extremized. The direct methods ■ 
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are not investigated in this thesis, but because of their importance, 
a brief reference to some of these methods is made. The most pop- 
ular direct methods are the gradient or steepest ascent methods 
developed independently by Kelley [2], and Bryson, Denham, Carroll and 
Mikami [3,^]. Many extensions to the basic method have been made and 
are discussed in the recent text by Sage [ 5 ] • Significantly different 
direct approaches are Bellman's dynamic programing [6,7], the conjugate 
gradient method discussed by Lasdon, Mitter, and Waren [8]; and the 
Optimal sweep method introduced by McReynolds and Bryson [ 9 ]. 

Indirect methods are those in which the control histories are not 
directly manipulated. Instead, the control problem is transformed into 
a two-point boundary value problem by deriving certain ordinary dif- 
ferential equations and boundary conditions which must be satisfied for 
mathematical optimality. The governing differential equations and 
boundary conditions are obtained from either the necessary conditions 
of the classical calculus of variations (see for example. Bliss [l] or 
Sage [5]), the Pontryagin maximum principle [10], or the theory of 
dynamic programing as discussed by Dreyfus [ 11 ] . The various succes- 
sive approximation techniques for solving the resulting nonlinear two- 
point boundary value problem are called indirect optimization methods. 
There are several generally applicable methods which have been developed 
•in recent years. A discussion of. these methods is deferred until 
Chapter II, where the general boundary value problem to be considered 
in this thesis is presented. 
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A common problem with existing indirect optimization methods is 
that initial approximations to either the solution or initial conditions 
of the boundary value problem are required. The success in solving 
the problem is sometimes extremely sensitive to the accuracy of the 
initial approximations. Convergence, when it occurs, is generally more 
rapid for indirect methods than for direct methods. Ideally, one seeks 
a method which converges rapidly to a solution from an arbitrary initial 
approximation. One of the major reasons for this investigation is to 
improve upon the performance of existing indirect methods with regard 
to this ideal characteristic. To this end, several new innovations and 
improvements to known techniques are combined in a unified approach. 

This includes the introduction of a power series integration method 
which exhibits several characteristics uniquely suited for determining 
numerical solutions of nonlinear two-point boundary value problems. 
Moreover, a new approach is presented for solving problems where the 
terminal boundary conditions are general functions of the final state 
and unknown final time. The computational algorithm is derived such 
that the differential equations to be integrated have improved numerical 
stability. Consequently, numerical difficulties due to ill conditioned 
matrices of boundary values can be avoided. 

The indirect optimization method presented here is applied to the 
solution of a minimum time, planar, Earth-Mars transfer problem for a 

constant low-thrust rocket. The problem was chosen because it has been 
used to test many other methods and thus, a direct comparison with 



previous results could be made. The optimization method is also applied 
to a similar problem where a minimum time rendezvous with- a moving target 
(the planet Mars) is required. 



CHAPTER II 


THE INDIRECT APPROACH TO TRAJECTORY OPTIMIZATION 

In order to apply an indirect trajectory optimization method, it 
is necessary to formulate the optimal control problem as a two-point 
boundary value problem. This is accomplished by deriving certain ordi- 
nary differential equations with accompanying boundary conditions which 
must be satisfied. In this chapter, necessary conditions from the cal- 
culus of variations are used to formulate the control problem as a two- 
point boundary value problem and previous numerical methods for solving 
such problems are discussed. 

The control problem which will be considered in this thesis can be 
stated: a system is described by a set of nonlinear ordinary differ- 

ential equations 

X ( t ) = f(x(t) ,u(t) ,t) (2.1) 

where x is an n vector and u is a q vector, q < n. The initial 
state is specified at some initial time t 



x 
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and the terminal state x(t f ) and time t f are given implicity by 
the following k vector of functions 

- *(x( t f), t f ) = 0 (2.2) 

It is necessary to determine the q vector of controls u(t) from some 

admissible set of controls so that the performance' index <j>(x(t f ) ,t f ) 

is minimized. The admissible set of controls will be taken to be the 

set of all piecewise continuous functions on the interval [t ,t ]'. 

An augmented functional is formed by adjoining to the given per- 

» 

formance index, using Lagrange multipliers, the contraints (2.1) and 
(2.2) to obtain 

t f ■- 

J = *(x(tf),t f ) + V T ip (x^t f ) ,t f ) + r X T (t)(f - x) dt (2.3) 

t 

o 

where v is a k vector of constant Lagrange multipliers, and X(t) is 
an n vector of time dependent Lagrange multipliers. The problem can now 
be viewed as one of seeking a minimum of the performance index J sub- • 
Ject to no additional constraints. A necessary condition for J to 
have an extremum is that the first variation of J vanish. 

It is helpful when considering variations of J to introduce the 
scalar function called the Hamiltonian l _ 

H(X,x,u,t) = X^ f(x,u,t) 
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such that by definition 

i = f(x,u,t) = (||) 

Using the Hamiltonian and an integration by parts, the performance index 
is written 

J = H^f) jt f) + yT 


t f A 


- X T (t)x(t) + J |h + x T xjdt 


t t 
o o 


The requirement that the first variation of J vanish along an optimal 
trajectory with final time not specified results in the following neces- 
sary conditions. 


1 - (f) T ' <*■"> 

Ht) T ' <*•» 


3H 

8u 


0 


(2.6) 
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( 2 . 8 ) 

(2.9) 

The derivation of these equations is adequately treated in several 
texts and the reader is referred in particular to Bliss [l] for classical 
problem formulations or to the recent book by Sage [5] for a treatment 
in the more modern control notation. Equations (2.4) and (2.5) repre- 
sent 2n simultaneous ordinary differential equations in the 2n + q- 
variables of the vectors X, x, and u. Equation (2.6) provides q con- 
ditions by which the q variables of the control vector u can be elim- 
inated from equations (2.4) and (2.5) so that 2n differential equations 
in 2n dependent variables are obtained. It is assumed that this is 
possible, since, in general, it is not always possible to solve explic- 
itly for each of the control variables via equation (2.6). 

The conditions (2.4) to (2.9) are merely necessary conditions for 
an extremum of J. A further necessary condition for J to be mini- 
mized is given by the non-negativity of the Weierstrass E-Function 



E = F(t ,x,X,U,X) - F(t,x,x,u,X) - ~{t ,x,x,u,A ) (X - x) 


9F • 

- g^(t,x,x,u,X)(U - u) > 0 
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where F = X (f(x,u,t) - x) and X and U are nonoptimal but 
permissible values for x(t) and u(t). This condition is normally- 
applied in conjunction with equation (2.6). For the example problems 
which will be presented, this condition is used to resolve an ambiguity 
in sign resulting from the application of equation (2.6) (see Appendix A). 

Equations (2.7), (2.8), and (2.9) yield n, k, and one algebraic 
equations, respectively, which must be satisfied at the final time. 

Since the Lagrange multipliers v are constants, v = 0. Adjoining 
these k trivial differential equations with those of equations (2.U) 
and (2.5) yields the following 2n + k system 



v * 0 / 


( 2 . 10 ) 


with the 2n + k + 1 boundary conditions needed in the case of unknown 
final time, given by 


X N • *0 


(t\ = 

+ V T ii 

V V 

3x 

t 

: f t 


♦( x (‘f) ^f) - 0 

H l + f| =0 


(n conditions) 

(n conditions) 

> 

(k conditions) 

(l condition) 


( 2 . 11 ) 
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To simplify reference to these equations, the system (2.10) is 
written as an N vector of nonlinear first order differential equations 



• 

z = 

f(z,t) 

(2.12) 

and the boundary conditions 

(2.11) 

are generalized to 


= 

z . 

Ol 

i 1 ^ 2 y « • • m 

(2.13) 

. h t ( z ( t f)’ t f) = 

0 

i = 1,2. . .(N - m + 1) 

(2.14) 


Thus, the solution of the optimal control problem is reduced to finding 
the solution of a system of nonlinear ordinary differential equations 
with two-point boundary conditions , the terminal boundary conditions in 
general being nonlinear functions of the terminal state and unknown 
terminal time. 

It should be noted that for many problems the terminal conditions 
(2.8) may not be very complex, and, as a consequence, the Lagrange 
multipliers v may be eliminated from equations (2.7) and (2.9) ana- 
lytically 'so that N + 1 terminal boundary conditions not involving 
v are obtained. This allows the deletion of the k trivial differ- 
ential equations v = 0 with the associated k terminal boundary 
conditions, and thus reduces the dimension of the problem considerably. 
This approach is used in the example problems which will be presented, 
although the computational algorithm' which is developed in succeeding 
chapters can be applied to the more difficult problem given by equa- 
tions (2.10) and (2.1l). 
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Since the system (2.12) is, in general, nonlinear, with two-point 
boundary conditions (2.13) and ( 2 . lU ) , solutions are not easily obtained 
The essential problem involved is to determine the missing initial con- 
ditions for the Lagrange multipliers so that at some later time equa- 
tions (2.l4) are satisfied. 

A numerical method for solving the more simplified version of the 
above boundary value problem with fixed final time was considered in 
19^9 by Hestenes [12], who also formulated the general optimal control 
problem given above [ 13 ] - Hestenes [l4] explained that his early work 
was not actively pursued due to a lack of interest in the problem. 

Breakwell [15], in 1959, published the general control problem 
formulation in the form given above and presented numerical results 
for a variety of problems . The problems were solved by repeated numer- 
ical integrations of the nonlinear differential equations with perturbed 
initial conditions and using an interpolation scheme for determining 
the initial conditions which would yield the desired terminal condition. 
A similar approach was used by Melbourne [l6], and Melbourne, Sauer, . 
and Richardson [17] for solving fixed time duration optimal payload 
trajectories for continuous low thrust orbit transfer maneuvers between 
the Earth and several other planets. These efforts are representative 
of some of the first attempts to obtain solutions by straightforward 
"brute force" tactics. These methods resulted in considerable frustra- 
tion and generally poor convergence or no convergence at all. Although 
some success was realized through such approaches , the general problems 
with convergence motivated the development of the direct methods. 
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The systematic algorithms of contemporary indirect optimization 
methods can be traced to the papers of Goodman and Lance [l8] in 1956 
and the work of Kalaba [19] in 1959, although neither of the papers 
was directly concerned with trajectory optimization. Goodman and Lance 

[18] discussed numerical solutions of systems of linear differential 
equations with two-point boundary conditions by the adjoint equations 
of Bliss [20]. They also proposed a method called complementary func- 
tions which utilizes the principle of superposition of particular and 
homogeneous (or complementary) solutions. In addition, they outlined 
an approach for solving nonlinear problems by relating initial and 
final boundary value perturbations of a nominal solution with a system 
of linear adjoint equations. Kalaba [19] developed the early ideas 

of- Hestenes [12] and produced a method conceptually different from those 
proposed by Goodman and Lance [ 18 ] . This method was called Quasiline- 
arization and required iterative solutions of a system of linear dif- 
ferential equations which were derived from a Taylor series expansion 
of the nonlinear equations. An initial solution approximation 
which satisfied initial and final boundary conditions was iteratively, 
improved by repeatedly solving the derived linear equations. Kalaba 

[19] gave a convergence proof and demonstrated the method for second 
order differential equations. Both of these approaches for solving 
nonlinear boundary value problems were restricted to fixed intervals 

of -the independent varilibXe^and^imple boundary conditions, and, there- 
fore, were not directly applicable to the general boundary value problem 
considered here. Extensions of the methods soon followed, however. 



and in this thesis those stemming from the ideas of Goodman and Lance 
[l8] are called Perturbation methods while those following Kalaba and 
Hestenes are called Quasilinearization methods. The Perturbation and 
Quasilinearization classifications will be more clearly distinguished 
in Chapter IV. 

The adjoint equation Perturbation approach of Goodman and Lance [l8] 
was extended to solve variable final time optimization problems by 
Jurovics and McIntyre [21], Jazwinski [22] developed the method further 
to allow for boundary conditions which are general functions of the prob- 
lem variables and time. A procedure for handling inequality constraints 
on state and control variables was also presented. Breakwell, Speyer, 
and Bryson [23] independently derived a method similar to Jazwinski's 
through considerations of the second variation of the calculus of 
variations. 

The alternate Perturbation approach of Goodman and Lance [l8], 
involving complementary functions, was also studied by Breakwell, Speyer, 
and Bryson [23] and compared to the adjoint Perturbation method from an 
operational standpoint of a computer storage requirement versus a matrix 
inversion requirement. Lewallen [2U] made extensive comparisions of the 
two Perturbation techniques and found them to have equivalent convergence 
characteristics. Further study of the Perturbation methods have been 
made by Shipman and Roberts [25] and Lastman [26] to show their connec- 
tion with the famous Kantorovich theorem [27] on Newton's method in 
functional analysis. Armstrong [28] has proposed a Perturbation method 
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which seeks to iteratively reduce a norm of terminal constraint dis- 
satisfaction and which displays some characteristics of direct methods. 

Adaptation of the Quasilinearization approach for optimal control 
problems was studied by McGill and Kenneth [29], [ 30 ] , who extended 
Kalaba's [193 convergence proof for systems of differential equations, 
and modified the method to solve variable final time problems. Their 
approach for solving variable final time problems involved the solution 
of a sequence of fixed final time problems and was inefficient. 

A novel approach for solving variable final time problems with the 
Quasilinearization method was developed independently by Conrad [ 31 1 and 
Long [32] and involves a change of independent variable to one integrated 
between fixed limits. Further extensions and improvement of the change 
of variable approach have been proposed by Johnson [33] and Leondes and 
Paine [3^]. Leondes and Paine [3*+] have also extended McGill and 
Kenneth's [29] convergence proof for problems with bounded control vari- 
ables. A different technique for handling variable final time problems 
with the Quasilinearization method has been proposed by Lewallen [ 35 ] - 
This approach is similar to the one used by Jazwinski [22] for the 
adjoint Perturbation method, and Lewallen [35] has shown this method to 
have convergence properties superior to the other above-mentioned Quasi- 
linerization methods. This method is also applicable to problems with 
general-type boundary conditions. Numerical techniques for handling 
inequality constraints on control and state variables with the Quasi- 
linearization method have been studied by Kenneth and Taylor [ 36 ] and 
McGill [37]. 
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Although the methods of indirect trajectory optimization are well 
developed, the methods sometimes are unable to converge to a solution 
from arbitrary initial solution guesses. Van Dine [38], [39] has sought 
to circumvent this problem by solving the linear boundary value problem 
of the Quasilinearization approach with a finite difference technique. 
Results have been obtained by this approach for fixed final time prob- 
lems and control variable inequality constraints, but it is doubtful 
whether the accuracy of other indirect methods can be obtained. Although 
the method is claimed to avoid the convergence problems of other in- 
direct methods, no direct comparisons on convergence have been published. 

The comparison of various direct and indirect trajectory optimi- 
zation methods by Kopp and McGill [4o], Moyer and Pinkham [4l], Tapley 
and Lewallen [42] and Tapley, Fowler, and Williamson [43] have pointed 
out the desirability for an indirect method with ability to converge 
from poor initial solution estimates . These studies have indicated 
that direct methods are more likely to converge from poor solution 
estimates , but that indirect methods have more rapid and accurate con- 
vergence when it occurs. Various strategies have been suggested for 
improving the range of convergence of indirect methods, but implemen- 
tation of these strategies often requires considerable skill and effort 
on the part of the user in order to retain the rapid convergence char- 
acteristics of the methods. Several of these schemes have been investi- 
gated by Lewallen, Tapley, and Williams [44]. In spite of notable 
improvement with these strategies, convergence sensitivities remain a 


problem. 
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In the following chapters , a method for solving the nonlinear 
boundary value problem is presented, which displays convergence prop- 
erties superior to previously published indirect optimization methods. 
Both Quasilinearization and Perturbation approaches are considered, and 
a Perturbation approach is selected because of its minimum storage re- 
quirement, ease of implementation, and fewer necessary integrations per 
iteration. The method, which is not developed according to standard 
perturbation formulations, reveals a new scheme for handling the variable 
final time problem resulting in a few number of iterations required 
for convergence. Numerical difficulties which sometimes occur with 
adjoint equations or perturbation equations are avoided through an alter- 
nate method for solving linear boundary value problems. A power series 
numerical integration scheme is used which allows for a variable inte- 
gration step size and simultaneous integration of reference and perturbed 
solutions. This eliminates the approximations of functions evaluated 
on the reference trajectory necessary without simultaneous integration 
of reference and perturbed solutions. The characteristically high 
accuracy capability of power series integration, together with elimina- 
tion of approximations used in the iterative solution process, give the 
method presented here a capability for obtaining extremely accurate 
numerical solutions of boundary value problems in ordinary differential 
equations. ~ ” - • 



CHAPTER III 


SOLUTION OF LINEAR DIFFERENTIAL EQUATIONS WITH NONLINEAR 
TWO-POINT BOUNDARY CONDITIONS 

An integral part of the method presented in Chapter IV for solving 
nonlinear boundary value problems requires 'numerical solutions of linear 
differential equations with nonlinear boundary conditions. Numerical 
methods for solvihg linear differential equations with linear boundary 
conditions are well known and include (l) the method of complementary 
solutions [l8], (2) the method of adjoint equations [l8] , [20], and, 

(3) the method of Green's functions [45]. An alternate method which has 
received attention in several recent papers [16] to [ 50 ] is known as the 
method of particular solutions. The method is extended here in order 
to solve systems of linear differential equations subject to two-point 

nonlinear boundary conditions with an unspecified terminal value of the 

i 

independent variable. 

The method of particular solutions is very similar to the method 
of complementary functions with the exception that the general solution 
is obtained by superposition of several particular solutions of the given 
set of differential equations rather than superposition of a single 
particular solution and several complementary or homogeneous solutions. 
When numerical solutions with digital computers are to be obtained, the 
method of particular solutions displays several important advantages over 
the above-mentioned methods. First of all, unlike the other methods, 
only one set of differential equations (the given set) need be programed 
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for solution which reduces the programing complexity. More important, 
however, is the fact that each solution integrated is a physically 
possible solution. Therefore, each equation integrated possesses the 
stability inherent in the physical system model with the result that 
solution values at the boundaries are closer in magnitude than would be 
expected for values of homogeneous solutions. This generally results 
in more numerical accuracy in the determination of superposition con- 
stants from inversion of matrices of boundary values. The first stated 
advantage motivated Miele's work [ 46 ]. Holloway [ 47 ] encountered numer- 
ical instabilities with the method of complementary functions and was 
led to study superposition of particular solutions because of the second 
stated advantage. 

Other discussions of the method and various applications to two- 
point and multipoint boundary value problems in ordinary and partial 
differential equations are given by Luckinbill and Childs [ 48 ], Baker 
and Childs [49], and Heideman [50]. These applications have been 
limited to problems with linear boundary conditions at specified bound- 
ary points. A more general approach for solving problems with terminal 
boundary conditions given as general nonlinear functions of the problem 
variables and an unspecified terminal time is developed below. 

Consider the N dimensional linear vector differential equation 

_ _ . y_(t)_= A(tly _+_ b.(.t)__ . _ (-3-.1-) — 


o 



20 


where A and b are a given N * N matrix and N vector, respectively,, 
of time varying functions. Boundary conditions are given at the initial 
specified time t in the form 



i = 1,2. . .m < N 


and terminal conditions are specified as general functions 


(3.2) 


^i^y^tf) ,tf) =0 i = 1,2. . .r < N 


(3.3) 


If the terminal time t f is specified, then r = N - m. If t^ is not 
specified, then r = N - m + 1. 

A general solution of equation (3.1) satisfying equations (3.2) 
and (3.3) can be represented by 

) 

S+i 

y(t) = 52 a v P k (t), S = N - m 
k=l 


with the auxiliary condition 


S+l 



where any S of the p are linearly independent particular solutions of 
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equation (3.1), and the ct^ are superposition constants. Initial 

1c 

conditions for the p (t) are chosen such that 


“iCO = y oi 

i = 1,2. . .m 

p^t^ = any value 

i = m+1 ,m+2. . .N 

4'(0 ■ 6 iA}c) + Vk 

k = 2,3. . .(S+l) 
i = 1,2. ..N 


where 


( 3.10 


6 ik = 


1 if i i- k+m-l 


8. if i = k+m-l 


r ik 


0 if i # k+m-l 
y^ if i = k+m-l 


and p 


W4 < 


The particular choice of 8 ik and y^ k given, insure the condi- 
tion for linear independence of particular solutions. The choice of the 
constants 8^ and y^ is free except that 8^, y^ # 0, and 8^ # 1. 
These constants can he chosen, depending pn the sensitivity of the system, 
to control the magnitude of terminal values of the particular solutions . 
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The superposition constants are to be chosen so that the 

equations in (3*3) are satisfied and the condition for superposition of 
particular solutions 

S+l 

E a w - 1 = 0 ‘ (3>5) 

k=l 


is satisfied. In order to determine the a , a formal substitution of 
S+i • 

E a vP (■&) is made for y(t) in equation (3.3), and equation (3.5) is 
k=l . . 

written as ^ r+ -^» to obtain an r+1 vector of functions h with elements 


h i( y (V a 2‘* , * 0l s+l’ t f)* t f) = 0 1 

h r+l( a l’“2' ' ‘ a S+l) = ° 


(3.6) 


For the case where t f is unknown, r = S + 1, and equation (3-6) 
represents S + 2 nonlinear equations in the S + 2 unknowns a and t . 

K X 

The equations in (3.6) can be solved for the and t^. by a Newton- 

Raphson [51] iterative procedure. 


% 
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The Newton-Raphson procedure is employed in the following manner. 

An initial guess for the a, and t^ is written as a column vector 

k f 


( 0 ) 


S+l 


Successive approximations for the proper values of this vector are ob- 
tained by repeated solution of the. following equation 


» (W1) = a (n) - [j (n) ] _1 h (n) 


n=0 ,1,2... 


(3.7) 


where is the Jacobian matrix with elements 






dh. 

3 

da , 


dh. 

l 

dt „ 


1.2. . .5+2 

1.2. . .5+l 


1,2. . .S+2 
S+2 


(3.8) 


and the n superscript denotes evaluation with the nth approximation for 

An) 



2b 


Expanding the derivatives of equation (3.8) and denoting elements 
of y(t) by y % 


fHi fulfil 

da j ' a 8 * t 8 “3 


3h. N 3h. . v 


+ 


3h. 

i 



since 


Also, 


y t (t) 


Y.- B.pjft) 

j=l J 


dh. N 9h. 3h. N 3h. 

"7 * £ + “7 = £ *7 




The method of solution requires the forward integration of the 
(S+l) particular solutions of equation (3.1), p (t), with initial con- 
ditions given by equation (3.U) from t Q to some assumed final time 
t . At the assumed final, time, the Jacobian matrix and the functions 
h. are evaluated using initial guesses for the a. . The equation (3.7) 

1 V 

yields new approximations for and t^. To continue the iteration 

of equation (3.7), it is necessary to integrate the particular solutions 
from the assumed final time to the new estimate for final time. The 
forward and backward integrations which may be necessary from the se- 
quence of final times generated by the iteration of equation (3.7) may 
be excessively cumbersome for same numerical integration methods. How- 
ever, the power series integration scheme discussed in Appendix B is 
well suited for this problem. 
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Application of the power series integration method yields Mth order 
polynomial approximations for p^(t^) and P^(tf) written as 



(3.10) 


where the a and b . are known power series coefficients 

A ) A> | 1 

determined by the method of Appendix B and t is used as the origin 

Q» 

of the power series expansions. If t is sufficiently close to the 

o - a 

true final time, the equations (3.10) represent sufficiently accurate 
• formulae in the application of equation (3.7). If | (t - t )| becomes 
too large, a new center of expansion can be used as explained in Ap- 
pendix B so that a specified accuracy is retained. 

With sufficiently close initial approximation for arid t^> 

the sequence (3.7) is rapidly convergent and yields the desired values 
for the and t f . Upon convergence of equation (3.7), the general 

solution of equation (3.1) satisfying equation (3.2) and equation (3.3) 
can be obtained by integrating (3.1) over the interval [t Q ,t f ] with 
initial conditions 


y.- 


i(* 0 ) = y oi 


i = 1,2. . .m 


S+l 


y i( l o) ' £ Vi ft) 


i = m+l,m+2, . . .N 
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In this manner, the solution y(t) can be constructed without storing 
the solutions p (t). 

Although the method requires initial approximations, this has not 

been found to be a problem. If it happens that estimates are available 

for the missing initial conditions, then using these estimates for 

p^t ) implies a, should be chosen near unity with other values of 
o 1 

near zero. A very accurate method of estimating the is obtained 

when the method is used in an iterative technique for solving nonlinear 
boundary value problems. This is discussed in the following chapter. 



CHAPTER IV 


METHODS FOR SOLVING NONLINEAR TWO-POINT BOUNDARY VALUE 
PROBLEMS WITH NONLINEAR BOUNDARY CONDITIONS 

A fundamental idea in most contemporary methods for solving non- 
linear boundary value problems in ordinary differential equations is to 
iteratively solve an associated set of linear differential equations. 

The solutions either converge to the nonlinear solution or provide a 
sequence of initial conditions which converge to the proper set of 
initial conditions for the nonlinear system. The linear differential 
equations are obtained by Taylor Series expansions of the functions 
defining the nonlinear system. This linearization process is common to 
many methods appearing under the various titles of Quasilineariza- 
tion [19], [31], [35]; Generalized Newton-Raphson Method [ 30 ], [ 32 ]; 
Modified Newton's Method [ 26 ]; Second Variation Methods [23]; Adjoint 
Method [21]; and Method of Perturbation Functions [42]. The systems 
of linear differential equations used by these methods are similar and 
in many cases identical. However, the actual sequence of approximate 
solutions generated may differ considerably depending on the type of 
initial solution approximation used, the manner in which given boundary 
conditions are employed, and the reference solution used in the linear- 
ization expansion for each iteration. 

The type of reference solution used for each iteration provides a 
basic classification of the methods into two groups which in this thesis 
will be called Perturbation methods and Quasilinearization methods. 
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Perturbation methods use a solution of the given nonlinear system with 
approximate values for unknown initial conditions as the reference solu- 
tion. The reference solution for Quasilinearization methods is a solu- 
tion of a system of linear differential equations derived from the 
nonlinear system. Both Quasilinearization and Perturbation methods are 
developed below using the ideas of Chapter III for obtaining solutions 
of linear systems with nonlinear boundary conditions. The Quasilineari- 
zation method obtained is recognized to be very similar to one proposed 
by Lewallen [35]. The Perturbation method obtained is significantly 
different from previously derived Perturbation methods and offers some 
decided advantages over presently known methods . 

Consider the system of N nonlinear differential equations 

z = f{z,t) (4.1) 

with initial conditions 

z . / 1 ] ■ z . i = 1,2, ...m (4.2) 

i\ o/oi 

and final conditions and final time given implicitly as the first time 
the following q = N - m + 1 vector of functions is satisfied 

bfMM’tf] = 0 ■ 

In general, a solution of equation (4.1) with initial conditions, equa- 
tion (4.2), and arbitrary values for the unspecified initial conditions 
will not satisfy the terminal conditions, equation (4.3). Denote such a 
solution by ^z(t), where the superscript is used to index the first 


0 
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approximation to a solution of equation (4.1) satisfying both the initial 
and final conditions. 

Let v(t) be an N vector of functions satisfying 
"l^o) * 2 oi i = 1,2. . .m 


h [ w ( t f)* t f] = 0 

as well as the differential equation 

w = f(w,t) t < t < t. (4.4) 

o ~ - f 


where the vector of functions f in equation (4.4) is the same vector 
of functions appearing in equation (4.1). Clearly, w(t) is a desired 
solution of equation (4.1) satisfying both equations (4.2) and (4.3). 

The functions appearing on the right-hand side of equation (4.4) can be 
expressed in a Taylor series expansion about the reference or approximate 
solution 1 z(t). That is 


w s f(w,t) 


(i .] sfrz.t) 

( 1 ) 

\ Z - t > * "Vz 

U - z / 


(4.5) 


If the expansion is truncated after the second term, the following 
linear differential equation, which is an approximation to equation (4.5), 
is obtained 


1 y = A(t) 1 y + b(t) 


w « 


(4.6) 
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where 



and 


b(t) = f( 1 z,t) - A(t) 1 z(t) 

A solution of equation (4.6) subject to the boundary conditions, 

= 2 oi i* 1 - 2 ---" 

»'[ 1 y( 1 t f ) = o (k.8) 

will yield an approximation for w(t). The solution ^y and correspond- 
ing value for final time H can be obtained by the method of particu- 
lar solutions described in Chapter III. The accuracy of the approximate 
solution obtained depends on the closeness of the solutions and w. 

Assume for the present that ^"z is sufficiently close to w so that 
V is a better approximation than 1 z. The manner in which additional 
approximate solutions are obtained determines whether the approach taken 
is categorized as a Perturbation method or a Quasilinearization method. 

A Quasilinearization Method 

Since ^y is a time varying vector of functions which is a better 
approximation for w than is ^z , then it is reasonable to replace ^z 
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with 1 y in the Taylor series expansion of equation (4.5). The linear 
differential equation (4.6) is then written 

v « 2 y w f^y.t) + ( 2 y -V) (4.9) 

and a solution of this equation subject to the given boundary conditions 

2 

provides a new approximate solution y. According to the definition 

set forth previously, this is a Quasilinearization method, the reference 

solution in the Taylor series expansion being a solution of the linear- 

2 

ized differential equation. Using y as the next reference solution 
yields ^y and so on for ^y, "V, ^y,... n y* Under appropriate condi- 

tions provided by the convergence proofs of references [29] and [34], 
the sequence of solutions n y converges to the desired solution w. 

These convergence proofs are restricted to problems with more simple 
boundary conditions than those expressed in equation (4.3). 

A closer observation of the Quasilinearization method with regard 
to the technique presented in Chapter III for ..solving the linear system 

with nonlinear boundary conditions reveals a fundamental difficulty. It 

2 2 
may happen that t^., the terminal time corresponding to the solution y, 

may be larger than the value of t obtained in the solution 

for '''y. In this case, \ is not known beyond "^t f , and the differen- 

tial equation (4.9) is not defined over the necessary range t < t < t f . 

However, if only one'Newtdn-Raphson iteration with equation (3.7) is made, 

and a linear extrapolation for V is used on the interval [VS-]- 

it is not necessary to integrate ^y past 1 t^ > in order to construct 
2 

y and a workable iterative scheme is realized. Lewallen [35], using 
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somewhat different arguments, has presented a Quasilinearization method 
which solves the variable final time problem essentially in the same 
manner proposed here. The primary difference between the method sug- 
gested here and the one proposed by Lewallen is that particular solutions 
rather than homogeneous solutions are used to construct the general solu- 
tion of the linear system. 

Lewallen [35] has compared this Quasilinearization approach with 
the Quasilinearization techniques described by Long [ 32 ], and McGill and 
Kenneth [30], and has found this particular approach to have better con- 
vergence characteristics. Another attractive, feature of the method is 
the capability for solving problems with general terminal boundary con- 
ditions of the form given in equation (4.3). 

A Perturbation Method 

In order to proceed with the development of the Perturbation method, 
it is assumed that ^y (the solution of the problem defined by equa- 
tions (4.6) to (4.8)) has been obtained. It is further assumed that ^y 
is a better approximate solution than ^z. (A method for satisfying 
this assumption is discussed later in the presentation.) In particular, 
^y(t Q ) should be a better approximation for w(t Q ) than was ^(t^). 

It is reasonable to expect that a solution of equation (4.l) using ini- 
tial conditions, ^y(t ), will be a better approximation for w than 

1 2 12 
was z. Such a solution is denoted by z. Replacing z with z 
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in the Taylor Series expansion (4.5) and the associated linear differen- 
tial equation (4.6) yields 


. ^ 2. (2 A 

3f( 2 z.t) 

(2 2 \ 

w « y ■ f^ z,t) 

3z 

( y - z ) 


(4.10) 


If the solution y is to closely approximate w, it must also satisfy 
similar boundary conditions , hence 




(4.11) 


/ 

2 

The solution y can be obtained by the method described in Chapter III, 

and the difficulty with final time encountered with the Quasilineariza- 

2 

tion method can be avoided. This is accomplished by generating z by 
integration of equation (4.1) simultaneously with the integration of 

all particular solutions of equation (4.10). In this manner the refer- 
2 

ence solution z is integrated to each estimate of final time required 

2 

by the algorithm of Chapter III. The solution z is thus defined over 
the entire range of time required for the solution of equation (4.10) sub- 
ject to the boundary conditions (4.1l). Consequently, as many Newton- 
Raphson iterations as desired with equation (3.7) can be made. This gives 
one the capability to obtain the initial conditions y(t Q ), so that the 

terminal conditions (4.11) are satisfied as accurately as desired. This 

capability has not been possible with existing methods for solving the 

2 

problem defined by equations (4.1) to (4.3). Once y(t Q ) determined, 

3 

a third reference solution z can be generated by integration of 
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equation (U.l) using initial conditions y(t o ). The new reference solu- 
tion allows continuation of the iterative process. If ^z is suffi-. 
ciently close to w, the sequence of initial conditions n y(t Q ) converges 
to w(t Q ), and hence, n z converges to w. By the definition set forth 
previously, this process is a Perturbation method, since the reference 
solution at each iteration is a solution of the given nonlinear system. 

Modifications for Improving Convergence 
In the foregoing discussions of Quasilinearization and Perturbation 
methods, it was assumed that the starting solution ^z was sufficiently 
close to the true solution w, so that convergence of the methods resulted. 
In practice, it is sometimes difficult to find an initial approximate 
solution which will lead to convergence. Consequently, various modi- 
fications of the basic methods outlined above have been proposed to 
improve convergence characteristics [23], [26], [3 1 *], and [1+4]. These 
modifications are commonly referred to as "iteration schemes" or "con- 
vergence schemes," and the usefulness of a given method is often closely 
tied to the "convergence schemes" employed. Two modifications to the 
basic procedures described above are discussed here. They should be 
considered to be integral parts of the basic methods rather than 
schemes which are Just added on. 

Consider a simple version of the nonlinear problem described by 
equations (4.1) to ( b.3 ) such that the terminal time is specified and 
terminal boundary conditions are given in the form 

fi 




i = k + l,k + 2,...(k + N - m) 
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for some k. (The restriction of fixed final time in this discussion 
can he removed by employing the transformation of independent variable 
as described by Conrad [ 31 ] and Long [32].) The following discussion 
assumes that the Perturbation method is used, but parallel arguments 
can be made for the Quasilinearization method. 

An approximate initial solution ^ z is obtained in the manner 
described previously to yield at the fixed final time the values 



Let w be a solution of the problem. Instead of, as before, seeking 
an approximation for w on each iteration, a function n w (in this case 
n = l) is sought which has initial conditions 

^ifo) = Z oi 1 = 1 

and terminal, conditions 



i = k + l,k + 2,...(k + N - m) 

0 < e < 1 


That is, w has terminal values between those of the reference, 
solution and the. desired solution. By choosing e sufficiently small, 
the approximation 


n* 

w 



n 


z 


) 


(i*. 12) 
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with 


ny i(tf) = ^i^f) i = k + l,k + 2,...(k + N - m) 

can be made as accurate as desired. To insure convergence, the approxi- 
mation must be sufficiently accurate such that the initial condition 
vector 


n+1 H‘o) ■ 

is sufficiently close to n w(t Q ), and hence, n+1 z(t f ) sufficiently 

close to n v(t ), so that n+1 z(t J is closer to v(t_) than is n z(t„) 
f f f f 

By choosing e sufficiently small, initial conditions for successive 
reference solutions are obtained which yield terminal conditions closer 
to the desired conditions than the previous reference solution. Repeti- 
tion of the process for n = 1,2,3... results in the construction of 
a sequence of terminal conditions n w(t^,) converging to w(t^), a 

sequence of initial conditions n z(t ) converging to n w(t ), and hence, 

o o 

solutions n z converging to w. 

A practical consideration is that some method for choosing z is 
necessary. From the above discussion it is apparent that there exists a 
value at each iteration which will work, but no a priori means of deter- 
mining e at each iteration is known. A trial and error method could 
be employed but this could be very inefficient. A simple method for 
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choosing e is proposed here which is based on practical considerations 
and which improves chances for convergence in a specified maximum number 
of iterations. Since it is apparent that the numerical methods proposed 
here will be used only with the aid of digital computers, it is always 
necessary to limit the number of iterations which can be attempted in 
order to conserve computer time and costs. 

First it will be noted that the scheme given above for choosing 
boundary values is equivalent to the following procedure. 

1. Solve equation (4.12) subject to the given initial and terminal 
boundary values 

z . i = 1 ,2, . . .m 
01 

n y^t^ = i = k + l,k + 2,. . . (k + N - m) 

2. Using the method of Chapter III, determine trial values for 
the missing initial conditions 

n y^t^ i = m+l,m+2,...N 

3. Compute a final set of missing initial conditions for the next 
reference trajectory n+ ^z from 




~ ■ n, i( v o)' ^fV( l o)~-‘ n2 i 



i = ra + 1 ,m + 2,...N 


n+ ^ 

It is shown in Appendix C that the values -obtained for z^(t Q ) us ^ n 8 
a given e are equivalent for both the case where terminal conditions 
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axe modified before solving for initial conditions and the case described 
in steps 1, 2, 3 above. Implementation of the second approach is easier 
than the first approach described because the change in initial conditions 



can be computed before choosing z. Furthermore, the latter method is 
simpler to apply for the more general problem with unspecified final 
time and nonlinear terminal boundary conditions, since it is only nec- 
essary to modify the missing initial conditions computed for each new 
reference trajectory. 

The following method for choosing z is proposed. Let M be the 
maximum number of solution iterations which can be allowed because of 
limitations of computer time and costs. When initial guess values are 
chosen for the missing initial values of the starting solution, also 
estimate the maximum deviation of a guessed value from its true value. 
Denote the maximum of the estimated deviations by d. Also choose a 
suitable norm for measuring the computed change y(t Q ) - n z(t o ). For 
example 




Compute a maximum allowable "initial condition change step size" from 
cl 3 d - 

6 « — , for example, 6 = On each reference trajectory iteration. 
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compute the norm, P [ n y(t Q ) , n z(t o ) ] , and compare . its value to 5. 
Choose e according to 


<5 

p [‘ y ( t =)' nz ( t =)] lf p [ Dy ( t °)' nz ( t o)] ’ 5 
r if p[ n y(t 0 ) ,\(t j] < t 


(4.13) 


Initial conditions for successive reference trajectories are computed 
from 



i ■ m + 1, . . .N 


(4. Ik) 


For variable fined time problems it may also be necessary to modify 
successive final time estimates according to 


n+1. 


n 


t. + e n At 
f f 


(4.15) 


where n At f is a computed change for final time determined in the 

solution for n y(t ). 

o 

The attractive features of this method for inducing convergence 
from poor initial estimates is that (l) very little effort on the part 
of the user is required (all he must do is choose 6), (2) as the solu- 
tions begin to converge, the scheme does not retard convergence and 
(3) chances for convergence with one computer run are - maximized con- 
sistent with available computer time. Of course, if 6 is chosen to 
be much smaller than necessary, the rate of convergence may be slowed 
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considerably. On problems which do not display convergence sensitiv- 
ities, 6 should be chosen arbitrarily large so that rate of conver- 
gence is not retarded. 

In situations where 6 is chosen to be too large, a second modifi- 
cation .can be employed which may induce convergence. When 6 is too 
large, a typical behavior is for initial conditions to be chosen in the 
proper direction for several iterations but then as the true values are 
approached, one or more of the unknown values may oscillate about its 
true value on successive iterations. When this occurs, halving the 
computed change for the oscillating value will bring it closer to its 
true value. Thus, the following procedure is proposed. 

n+ 2 \ 

(a) Compute z i(t 0 ) as descr iked above, | 


i=m+l,m+2...N 

(b) Compute [“*\(t 0 ) - %( t „)]/[\( t 0 ) - n "\(\,)] 

(c) If the above quantity is less than -1/2, compute 


>( 4 . 16 ) 


n+l 


N - \( l o) - ItNw - °’\n] 


If the quantity in (b) is positive, the particular element of the vector 

is not oscillating on successive iterations. If the quantity in (b) is 

negative but larger theui -1/2, then the oscillation has a convergent 

nature. In either case there is no reason to modify the computed value 
n+l 


for 


z^(t Q ). This modification may also be applied to successive 


final time estimates for variable final time problems. 

There are other possible variations of the two basic modifications 


presented above. For example, one might reduce the value of 6 whenever 
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the norm p[ n y(t ), n z(t )] is computed to be less than 6, and/or 6 
o o 

might be reduced whenever oscillation of one or more initial condition 
values or final time value occurs. Details of such procedures are best 
worked out through numerical experiments. When upper and lower bounds 
are known for missing initial conditions and/or final time, these bounds 
should be imposed in the event that the values chosen violate these bounds . 

Comparison of Quasilinearization and Perturbation 
Computational Requirements 

A basic goal of this investigation is to formulate an improved com- 
putational method for solving nonlinear two-point boundary value problems. 
While convergence characteristics are a major concern, other factors 
such as ease of implementation, computer storage requirements, computer 
time per iteration, and control of solution accuracy are also important. 
Two basic methods. Quasilinearization and Perturbation, have been pro- 
posed from a theoretical standpoint. A comparison of the computational 
requirements and restrictions of each method is made here. This com- 
parison reveals the Perturbation method to be a more efficient computa- 
tional scheme, especially when used in a unified approach with the 
particular solution method of Chapter III and the power series numerical 
integration method discussed in Appendix B. 

A distinctive computational feature of the Quasilinearization 
method, often considered- to-be an advantage of 'the' method, 'is that if 
is not necessary to program the given nonlinear system of equations for 
solution. For convenience in the previous presentation of the method. 
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it was assumed that an initial guess solution was obtained by integra- 
tion of the given nonlinear system with assumed initial conditions. 

This is not necessary since any guess solution satisfying only the bound- 
ary conditions can be used. This "advantage" of the method is lost, 
however, if a starting solution is generated by integration of the non- 
linear system. With the Quasilinearization method, one has an option of 
storing each particular solution of the linear system at each integra- 
tion step and forming the reference solution by the properly weighted 
sum of these solutions, or one may avoid the storage problem by. inte- 
grating the linear system with the proper initial conditions to form the 
reference solution. With reference to equation ( 1 +. 9 ), the latter approach 
still requires that the values for f( n y,t) and 8f( n y,t)/8z n y be 
stored at each numerical integration step. To simplify access to these 
stored quantities, one is forced to use numerical integration schemes 
which use a fixed integration step size. The choice of this step size 
is influenced not only by truncation error of the numerical integration 
scheme, but also by the required spacing of the stored quantities in 
order to achieve the necessary accuracy for the approximation of these 
functions along the reference trajectory. Thus, selection of integra- 
tion step size in order to achieve a specified final solution accuracy 
is not a routine matter. The restriction to fixed numerical integration 
step size could be a serious handicap for problems where considerable 
integration speed and accuracy are realized through frequent changes in 
integration step size. Many of the boundary value problems arising in 
optimal control theory (for example, those in interplanetary navigation) 
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have this property. The Quasilinearization method, with the minimum 
storage option, requires N - m + 2 integrations of the linear system 
(eq. (4.9)) at each reference solution iteration since N - m + 1 
integrations are required to determine the proper initial conditions, 
and then these initial conditions must be used in one additional inte- 
gration of the system to generate the reference trajectory. 

In comparison, the Perturbation method offers some unique computa- 
tional advantages. Since it is never necessary to generate the entire 

solution of the linear system (4.12), but only the initial conditions 

n y(t Q ), there is no need for storing perturbed particular solutions of 

the linear system. Furthermore, since the reference solution n z can 

be generated by simultaneously integrating equation (4.1) forward with 

all particular solutions of equation (4.12), the quantities f( n z,t) 

and 8f( n z,t)/3z n z appearing in equation (4.12) need not be saved. 

n 

They are merely computed from z at each integration step, used in 
all integrations of equation (4.12) for the integration step, and then 
discarded. With this procedure, variable step integration schemes may 
be used since there is no need to restrict end points of numerical inte- 
gration steps to coincide with previously stored information. 

The combination of simultaneous and variable step integration of. 
the nonlinear and linearized differential equations which is possible 
with the Perturbation method provides an additional advantage for this 
approach. The variable step capability allows one to use integration 
schemes which automatically determine an integration step size to yield 



a specified solution accuracy. The simultaneous integration of the non- 
linear and linearized equations not only eliminates storage of the func- 
tions f( n z,t) and 2 n z, but it also eliminates the necessity 

u Z 

for interpolation schemes used to convert discreet values of these func- 
tions into more accurate approximations over the integration step. Si- 
multaneous integration automatically provides the interpolation for 
these functions through the mechanics of the particular integration 
scheme used. With the variable step power series integration method 

discussed in Appendix B, Taylor series expansions of these functions are 

/ 

generated which yield an approximation accuracy equal to the desired 
integration accuracy. The automatic step size selection of this method 
also relieves the user of the burden of determining beforehand an 

acceptable integration step size. 

i 

In addition to the above-mentioned computational advantages of the 
Perturbation method over the Quasilinearization method, the Perturbation 
method requires one less numerical integration per iteration of a com- 
parable set of differential equations. This may not be immediately 
obvious since it has been previously indicated that N - m + 1 integra- 
tions are required to solve the linear system (4.12) and one integration 
of equation (4.1) is necessary to construct the reference trajectory. 
This totals to N - m + 2 integrations per iteration, but only 
N - m + 1 are required if the following observation is made. 

Theorem 1: A solution n z of the nonlinear system (4.1) is iden- 

tical to a particular solution of the linear system (4.12) if 
initial conditions of the two solutions are identical. 



Proof: Let p be a particular solution of equation (4.12) and let 

n z be a solution of equation (4.1). Let x be defined 

x{t) = p(t) - n z(t) 

which implies 

x(t) = p(t) - n z(t) 

Since p satisfies equation (4.12) and n z satisfies equation (4.1), 



Now this is a homogenous linear differential equation, and by 
hypothesis 

x ( t o) = KO - M*o) '■ 0 

For these initial conditions, it is well known (see, for example, 
Petrovski [52]) that the solution for x(t) is 

_ - _ . — - x(t)- = 0 - - - - - - - - 



which implies 
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p(t) = n z(t) 

and thus the proof is complete. 

Using this theorem, one of the N - m + 1 integrations of equation (H.12) 
can be eliminated since the reference trajectory n z can be used in its 
place. 

A further point of comparison of the computational requirements 

for Perturbation and Quasilinearization methods is concerned with the 

manner in which convergence is detected for the methods. For the 

Perturbation method, a direct indication of convergence is given when 

the reference trajectory satisfies the terminal boundary conditions to 

some specified accuracy, or when the change in the initial condition 

vector is less than a specified accuracy. However, with the 

Quasilinearization approach, the reference trajectory does not satisfy 

the nonlinear system until convergence has occurred. To determine when 

successive trajectory iterations are converging, it is necessary to 

n n 1 

compute some suitable norm p[ y(t), y ( t ) ] . The computation of this 
norm requires a comparison of the successive reference trajectories at 
each integration step and consequently requires additional programing 
and computer time. 

This comparison of the computational requirements and restrictions 
of the Quasilinearization and Perturbation methods indicates that the 
Perturbation method is somewhat easier to implement, and is better 
suited for adaptation with the method of particular solutions described 


) 
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in Chapter III and the power series integration scheme presented in 
Appendix B. Outlined below is a computational algorithm which combines 
these various concepts together with the proposed modifications for 
extending the range of convergence in a unified method for solving the 
nonlinear two-point boundary value problem in ordinary differential 
equations. 

The Particular Solution Perturbation Method 
To obtain an efficient computational algorithm utilizing the con- 
cepts set forth in this chapter and the preceding chapter, a study of 
the manner in which these various ideas are incorporated into an inte- 
grated framework is in order. Because the Perturbation concept is 
employed with the method of particular solutions, the algorithm described 
below is referred to as the Particular Solution Perturbation Method (PSPM) 

On each solution iteration, the PSPM requires a simultaneous for- 
ward integration of the given N dimensional nonlinear system 

n z = f( n z,t) (i+.lT) 

together with S forward integrations of the derived linear system 


n y n r ♦ f(°z,t) - -£ "*■*) n , 


where S = N - m and m is the number of specified' initial conditions 



i = 1,2,. . .m 
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The terminal conditions specified for the nonlinear system are assumed 
to be of the form 

h i[ Z ( t f)’ t f] = ° i = 1,2, ... (S + l) 
so that the linear system is to satisfy boundary conditions given by 


\{\) * Z oi i = 

h 1 [ n y( n t f ), n t f ] =0 i = 1,2,. ..(s + 1) 

Using the method of particular solutions, n y is expressed 

' . S+l 

n y(t) = I] a k p (t) 
k=l k n 

Subject to 


■ ( U . 19 ) 


l 


S+l 

a k = 1 

where any S of the n p k . are linearly independent particular solutions 
of equation (1+.18), and the are superposition constants. Theorem 1 


£ 

k=l 


n p 1 (t) = n z(t) 


is used to write 
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and the systems (4.17) and (4.18) are written 



(4.20) 


with m initial conditions for each solution provided by 



i = 1,2,. . .m 
k = 1,2,. .. (S + 1) 


and other boundary conditions 



to be satisfied by selection of proper values for and n t^. 

Since only m initial conditions for the solution p^ are speci- 
fied, the remaining S missing initial conditions are taken to be the 
best available estimates for these values. For n = 1 , the missing 
initial conditions are actually estimates , but for n "=“273 j4‘, , the 
initial conditions are provided by the algorithm in the manner described 
previously for n z(t Q ) (eqs. (4.l4) and (4.l6)). Initial conditions for 
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k 1 

the n P » k = 2,3»...(S + l) are determined from ^p according to 

the scheme of equation (3.4), 


k 

n P i 





+ Y 


ik 


i = 1,2,. . .N 
k = 2,3,. . . (S + 1) 


where 


6 ik = 


if i f k + m - 1 


[0^ if i=k+m-l 


} (4.22) 


Y ik = 


(0 if i^k+m-1 


(Y i if i = k + m - 1 and 


n#o)| 


and 0^ and are perturbation factors prescribed by the user in 

order to control the magnitude of deviations between the various partic- 
ular solutions. 

At each iteration of the PSPM, S + 1 vector differential equa- 
tions (U. 20) are integrated from t to the best estimate for t„, and 

o I 

the Newton-Raphson algorithm, equation (3.7), is used to determine 
values of and t^, which satisfy the boundary conditions (4.21). 

However, in order to efficiently incorporate this algorithm into the 
PSPM, the following observations are' made. Each iteration of the 
Newton-Raphson algorithm yields estimates of the superposition constants 
from which an estimate of the initial conditions 


S+l ,, 

n >to) * £ “k rf (‘o) 
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can be made. This estimate can be used to compute an estimate for the 
change in the initial conditions of the nonlinear system, n Az(t Q ), 
where 





S+l 





Using equation (4.22) and simplifying, the estimated change in initial 
conditions can be expressed as a function of the and perturbation 
factors 


SW - 0 1 ■ i.. 2 ' 

” lz i( t o) = °k[n p i( t o)( 6 l ' U + T ik] 


.m 


i = m + l,m + 2,...N 
k = i - m + 1 


A suitable norm for this estimated change p[ n Az(t Q )] can be computed 
and compeared to the maximum allowable norm for this change (the value 6 
appearing in eq. (4.13)). If 


•r-wi 


> 6 


then additional iterations of the Newton-Raphson algorithm are not use- 
ful since this would only serve to compute n Az(t Q ) to greater accuracy, 
with the subsequent application of the convergence modification (4. l4) 
wasting this effort. Therefore, in this situation, only one Newton- 
Raphson iteration should be made. When the norm of n Az(t Q ) is less 
than 6, then an indication that the PSPM is in the terminal stages of 
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convergence is obtained, and continued' iterations of the Nevton-Raphson 
algorithm can be expected to yield better estimates of the unknown 
initial conditions and final time. The effect of the number of Newton- 
Raphson iterations allowed for the case when p[ n Az(t Q )] is less than 
6 is a subject of investigation in the following chapter. 

When the final Newton-Raphson iteration is made on each reference 
solution, and the subsequent estimate of n Az(t^) is obtained, the 
modification (U.l6) is applied to yield a final value for n Az(t Q ). 
Initial conditions for the next reference solution are then computed 
from 


an'tyo) = n pl ( t o) + "“(‘o) 

In this manner, if convergence occurs, the initial conditions 
n P^(t Q ') converge to the proper initial conditions of the desired non- 
linear solution. Since n y(t) also converges to the desired nonlinear 
solution, the following result is obtained at convergence 



n 

z 





> 


This condition is satisfied if = 1 and = ... = a s+1 = 

Besides offering a simple and positive test for convergence of the PSPM, 
the above mentioned final converged values for the provide reason- 

able estimates for these values which are required by the Newton-Raphson 
algorithm. These estimates become increasingly more accurate as the 


PSPM converges. 
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In the next chapter, the convergence characteristics of the PSPM 
are investigated and compared with published results for other Perturba^ 
tion and Quasilinearization methods. The effects on convergence by the 
various modifications are investigated separately in order to evaluate 


the effectiveness of each. 



CHAPTER V 


DISCUSSION OF RESULTS 

In this chapter, the convergence characteristics of the Particular 
Solution Perturbation Method (PSPM) are investigated on two typical 
nonlinear boundary value problems which result from the formulation of 
an optimal control problem for solution by an indirect method. The 
problems are formulated from the same basic optimal control problem and 
differ only in the boundary conditions which are imposed. The basic 
control problem considered is the determination of the thrust vector 
control for a minimum time, planar, Earth-Mars, orbit transfer for a 
spacecraft with a continuously firing, low-thrust rocket engine. This 
problem was selected because it has been used to test several other 
optimization methods, and consequently considerable data were available 
from which a direct comparison of results could be made. 

The equations of motion for the thrusting rocket are formulated in 
heliocentric, polar coordinates where only the gravitational attraction 
of the Sun is considered (Fig. 1). In addition, it is assumed that the 
thrust vector of the rocket can be turned continuously and effortlessly 
so that the spacecraft is idealized as a point mass with negligible 
rotational dynamics. The nonlinear ordinary differential equations to 
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Figure 1.- Coordinate system. 


be solved for the determination of minimum time transfer trajectories 
are derived in Appendix A and include the spacecraft equations of motion 



Z,- = m = -c = f 
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and the associated Euler-Lagrange differential equations 



where T is the constant thrust of the rocket, c is the constant 
mass flow rate of the exhaust, GM is the gravitational constant of 

— 2 I — ^ 5 

A^ + Ag and cos g = -Ag/^ *1 + *2 ‘ 

Example Problem 1 

For the first example problem considered, it is required only that 
the spacecraft reach an assumed circular Mars orbit with zero radial 
velocity and tangential velocity equal to that of Mars. The final 
angle 0 is not specified. The known initial conditions are the posi- 
tion, velocity, and mass of the spacecraft as it leaves an assumed 
circular Earth orbit; the normalized value of one Lagrange multiplier; 
and a known zero value for the constant A^, which results from not 
specifying a value for ©(t^); 

hN - to) - 0 
M‘o) - to) - 1 
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MM “ to) - 1 

O 

MM - e (M - 0 
MM ■ «(M ■ 1 
MM “MM '- 1 
MM = H ' 0 ■ 

The normalization of the Lagrange multipliers and other system param- 
eters is discussed in Appendix A. The terminal boundary conditions at 
the unknown final time are 

h i[ z (M M - MM - 0 • 0 

h 2 [ Z (t f )’tf] = Z 2^t) “ °- 8l012 728 = 0 

h 3 [ z (M’M * MM - 5236790 ■ ° 

For this problem, the dimension of the vector Z is N = 9, with m = 7 
specified initial conditions and S+l=N-m+l=3 terminal condi- 
tions given since final time is unknown. The unspecified initial condi- 
tions are 


MM * MM 
MM = MM 


/ 
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■ Example Problem 2 

For the second example problem, it is required that the final 
spacecraft central angle e(t^) be equal to the central angle of Mars 
at the time of rendezvous. The central angle of Mars at the end of the 
transfer trajectory is computed from a known central angle of the planet 
at the beginning of the transfer, the constant angular velocity of the 
Mars about the Sun, and the time of flight, 

e (M ■ e M (g + J 

Thus, for this problem an additional terminal boundary condition is 
added to the set given for example problem 1, 

h ( 2 (‘f)' i f] = sW • '4*o) ■ = 0 

Since, in this case, the terminal value of 8 is constrained, X^ 

cannot be determined to be zero, and the initial and constant value for 

X^ is unknown. Therefore, for this example problem there are three 

unspecified initial conditions: X, (t ), X 0 (t ), and X, (t ). 

1 o 2 o 4 o 

Numerical Results for Example Problem 1 
The solution of example problem 1 provided correct initial multi- 
plier values and final time as follows: 


X l( t o) = -°- U ^865 
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x 2 (t o j = ^i. 07855 

t f = 3.319^37 (l time. unit = 53.132355 days) 

\ 

These values were obtained, using a relative error "bound of 10 - ^ with. 

the power series integration scheme of Appendix B. Convergence of the 

PSPM was detected "by requiring that the sum of the absolute values of 

-k 

the superposition constants and be less than 1 * 10 . For 

this problem, this convergence criterion was more demanding than requir- 
ing that the initial condition and final time changes be less than the 
specified convergence tolerance, since it was observed that changes in 
these values were about one order of magnitude smaller than values of 
a 2 and a^. All computations were made in single precision arithmetic 
(eight significant figures) on the Uni vac 1108 digital computer. Each 
iteration of the PSPM required approximately 2 seconds of computer time. 

In order to evaluate the convergence sensitivity of the PSPM to 
initial guess values for X^, X,,, and t^., the problem was solved many 

times using starting guesses which deviated from the true values by 
known percentages. The deviations from the true values were chosen in 
a systematic manner so that the data could be presented in the form of 
convergence envelopes. The convergence envelope shown in Figure 2 was 
constructed from all initial guess data having a final time error of 
-20 percent (a guessed final time less than the actual final time). 

The convergence envelope was formed by locating the percentage devia- 
tions used for the initial guess- values of the two Lagrange multipliers 
on a Cartesian coordinate grid. Each problem attempted was located on 



hundreds of percent 
error 


6o 



Figure 2.- Convergence envelope for -20 percent final time error. 
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the grid, by a small circle. Darkened circles represent initial guess 
values which did not lead to convergence in 20 iterations of the PSFM. 

Open circles containing numbers represent initial guess values for the 
multipliers which did lead to convergence, and the numbers in the cir- 
cles represent the number of iterations required. On divergent trials, 
typical behavior of the method was to successively select multiplier 
changes in the wrong direction on each iteration. Also shown in Fig- 
ure 2 is the boundary of a convergence envelope obtained for this prob- 
lem by Lewallen [2U], who investigated and compared several trajectory 
optimization methods. In reference [ 2U ] , similar sized convergence 
boundaries were presented for three methods; the Method of Adjoint 
Functions studied by Jazwinski [22]; the Method of Perturbation Func- 
tions discussed by Breakwell, Speyer, and Bryson [23]; and Lewallen's [35] 
Modified Quasilinearization Method. These methods typically required 
11 to 20 iterations for initial multiplier errors along the outer edge 
of the convergence boundary shown. The superior convergence characteris- 
tics of the PSPM are evident. 

Presented in Figures 3 and k are similar convergence envelopes for 
cases with 0 percent and +20 percent deviations in initial guesses for 
'final time, respectively. Also shown in these figures are typical con- 
vergence envelopes presented in references [ 2h ] and [1+2] on the same 

problem with the. thre.e .methods, mentioned previously.. _The_ superior con- 

vergence characteristics of the PSPM are again indicated by these data. 



hundreds of percent 
error 






64 


The most probable reasons for this marked difference in convergence 
characteristics of the PSFM and the other methods (which are quite simi- 
lar methods) are discussed below. 

For the data presented in Figures 2, 3, and 4, the value used for 
the maximum allowable initial condition step size norm, 6 of equation 
(4.13), was 0.5. Each time that the requested step size norm was less 
than 6, the value of 5 was set equal to the norm of the requested 
change. For significant errors in the initial values of the Lagrange 
multipliers, typical values for p[ n AZ(t o )], the norm of the requested 
change in ^(t o ), and t^ varied between 8 and 5000. This 

means that for some cases, the value of t used in the convergence 
modification of equation (4.14) was on the order of 1 * 10 ^ and the 
requested changes in A 1 (t Q ), \^(t Q ) , and t f were reduced by this 

fraction. In comparison, the various methods studied in references [24] 
and [42] were implemented with a fractional correction scheme which had 
the essential effect of halving the computed initial condition changes 
and final time changes on the first few iterations. With this scheme, 
for multiplier errors below the indicated boundary in Figures 2, 3, 
and 4, the first iteration yielded multipliers in the upper part of the 
envelopes. The PSPM also diverges in these upper regions due to sub- 
sequent multiplier changes being selected in the "wrong direction." 
However, for multiplier guesses in the lower half of the envelopes, the 
fractional correction computed for the PSPM was sufficiently small to 
prevent "stepping over" the solution. Had the PSPM been implemented 
with the fractional correction scheme of references [24] and [42], the 
PSPM convergence characteristics would have been similar to the 
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convergence characteristics of the methods studied in these references. 
However, not all of the desirable convergence characteristics of the 
PSFM can he attributed to this one item alone. It is expected that many 
cases which converged with the PSPM would not have if the convergence 
modification (4.l6) had not been used. 

A vivid illustration of the importance of the modification (U.l6) 
is presented by the data of Table I. These data represent the values 
of the Lagrange multipliers and final time estimates obtained on suc- 
cessive iterations with and without the convergence modification (^.16). 
The initial guesses correspond to multiplier errors of 0 percent and 
-50 percent with a terminal time error of -20 percent. The PSPM will 
never converge from these initial guesses without the modification, and 
with it convergence is obtained in eight iterations . Of the eight iter- 
ations, only three required application of the modification as indicated 
by the (H) symbol in the table for values affected by the halving 
feature. 

Another feature of the PSPM which contributed in part to the good 
convergence performance shown in Figures 2, 3, and 4 was the use of 
upper and lower bounds for t^,. Upper and lower bounds of H.t and 2.2 
were specified, and although these bounds were rarely approached, they 
were imposed in several instances. Since bounds on the Lagrange multi- 
pliers X^(t Q ) and ^2^0 were no ^ easily determined, no upper and 
lower bounds for their values were specified in this study. 

For those starting guesses indicated in Figures 2, 3, and L which 
did not lead to convergence of the PSPM, the typical behavior of the 





COMPARISON OF PSPM WITH AND WITHOUT CONVERGENCE MODIFICATION (4.16) 
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(H) indicates value was obtained from the halving feature of the modification. 
1.0E-1 indicates 1.0 * 10 \ 
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PSPM was to select both multiplier initial value changes to be in the 
wrong direction on each iteration. Usually after 20 iterations the 
magnitude of the initial values for A^(t Q ) and A^t^ were so l ar f5 e > 
that the effects of A^ on the solution of the Euler-Lagrange equations 
were negligible. Consequently, although increasingly larger values were 
obtained for A^(t Q ) and ^ 2 ^ 0 ^* ^eir ra ^i° remained almost constant 
and each successive reference trajectory was an essential repeat of the 
previous reference trajectory. With this type of behavior, it was appar- 
ent that the PSPM would not converge in any number of iterations for the 
particular choice of initial Lagrange multipliers. Any initial guesses 
for A^(t Q ) and Ig^o^ which were large in magnitude compared to 

A_(t ) = -1, caused the PSPM to generate very similar reference trajec- 
3 o 

tories on the first several iterations. However, in most cases, the 

multiplier changes on these first few iterations were made in the proper 

direction, and convergence resulted. This behavior suggests that the 

convergence space of the PSPM is boundless in the lower quadrants of the 

envelopes of Figures 2, 3, and 4 when a value of 6 is used which will 

prevent the method from "stepping over" the solution and selecting 

values in the upper quadrants of the convergence envelopes. 

The operation of the PSPM is illustrated graphically in Figure 5. 

Each arrow represents the change in the values of A., (t ) and A 0 (t ) 

1 o 2 o 

taken on each iteration.. Also, presented in tabulated- form -is the value - 
of t f at each iteration, the requested step size norm, the fractional 
reduction factor e, and the value of the terminal constraint norm 
obtained with the reference solution of each iteration. The initial 
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Figure 5.- Successive iterates of the PSPM. 


69 


value of 6 was taken to be 0.5, and. 6 was set equal to the norm of 
the requested change whenever this norm was less than 6. Figure 5 
illustrates that the proper direction for the multiplier change is 
chosen at each iteration. Similar results can be expected for much 
larger errors in the lower left quadrant of the convergence envelope. 

It is interesting to note the behavior of the terminal constraint norm 
of the reference solution at each iteration for the problem presented 
in Figure 5. Although the PSPM takes each step in the proper direction, 
the successive values of the terminal constraint norm initially decrease 
very slightly and actually increase for the fourth, fifth, and sixth 
iterations. This behavior suggests that fractional correction proce- 
dures utilizing the value of a terminal constraint norm of the reference 
trajectory may not work well on this problem. 

An Improved Method for Choosing 6 
Perhaps the most appropriate criticism of the PSPM as presented is 
the necessity for choosing a value for 6, the maximum initial condition 
change norm. If 6 is chosen too large, the method may "step over" the 
solution into the divergent region. On the other hand, if 6 is chosen 
too small, the convergence of the method is unduly retarded. For exam- 
ple, if 6 = 0.25 had been selected for the problem presented in Fig- 
ure 5, it would have taken 14 iterations to arrive at the same 
multiplier values obtained in seven iterations when <5 was chosen to 
be 0.5. In order to illustrate the sensitivity (or insensitivity) of 
the PSPM to the value chosen for 6, the problem of Figure 5 was solved 
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using 6 values of 0.2 5, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75, and 2.0. 

The results of this investigation are presented in Table II. 

TABLE II 

NUMBER OF ITERATIONS REQUIRED FOR VARIOUS VALUES OF 6 


6 

No. iterations 
required 

6 

No. iterations 
required 

0.25 

20 

1.25 

13 

0.5 

12 

i.5 

Diverge 

0.75 

14 

1.75 

10 

1.0 

19 

2.0 

Diverge 


These results indicate that for 6 = 1.5, the allowable change in ini- 
tial conditions was large enough to allow the method to "step over" the 
solution. The convergence with 6 = 1.75 was coincidental since the 
second iteration produced the same multipliers as the 7th iteration of 
the case when 6 = 0.5. 

The behavior of the PSPM shown in Figure 5 suggests an approach 
for making the selection of 6 a self-adapting feature of the method. 
When each successive initial condition change vector is taken in the 
same direction as the previous change vector, an indication that 6 
can be increased is obtained. This behavior can be detected by forming 
the dot product of successive initial condition (and final time) change 
vectors and computing the cosine of the angle between successive change 
vectors. When this angle is near zero, successive Lagrange multiplier 
and final time values lie very near a line connecting the initially 
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assumed values for these parameters and values of these parameters near 

I 

the true values (as indicated in Fig. 5)* By referring to the table of 
Figure 5, one may plot the requested step size norm as a function of 
distance moved along this "line" which we will designate as the "con- 
vergence path." The data of the table in Figure 5 are plotted in this 
manner in Figure 6. At convergence, the requested change is zero. An 
estimate of the distance to move along the convergence path in order to 
obtain multipliers and final time which yield a zero change (converged 
values) is obtained by estimating the point of intersection of the 
curve and the horizontal axis in Figure 6— The slope of this curve can 
be computed numerically by evaluating the successive requested norms 
and keeping track of the distance moved along the convergence path on 
successive iterations. Graphically, the estimated distance to move 
along the convergence path using the self-adapting approach is shown by 
the intersection of the dashed lines and the horizontal axis in Fig- 
ure 6. To implement this self-adapting approach, the initial iteration 
is made with any small value for 6, (6^ = 0.5 in Fig. 6). If the 
change vector of the second iteration is in the approximate same direc- 
tion as the first iteration, then the distance to move along the conver- 
gence path is found by 

(Requested norm 2)|6^J 

^2 (Requested norm 1) - (Requested norm '2 )' -- - 



Requested step size norm 
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- Behavior of requested step Bize norm along convergence path. 
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Successive values for 6 at each iteration are chosen in this manner 
until successive change vectors are not approximately along the same 
line or until the computed value for 6 is less than the initial spec- 
ified value <5-. 

This self-adapting feature was incorporated into the PSFM and 
studied on a variety of initial multiplier guesses. Typical results 
with this scheme are illustrated by solving the problem of Figure 5 
with <5^ = 0.25 in 12 iterations instead of the 20 required when 
5 = 0.25 at each iteration. Using <5^ = 0.25 with initial final time 
error of -20 percent and initial Lagrange multiplier errors of 
-1000 percent for both and convergence was obtained in 16 

iterations. Similarly, with initial multiplier errors of +1000 percent, 
-1000 percent, and final time error of -20 percent, convergence was 
obtained in 14 iterations . 

Experience with this scheme is limited at the present time and 1 
undoubtedly its effectiveness is somewhat problem dependent. For exam- 
ple, if the curve of Figure 6 were concave instead of convex, the 
scheme may cause 6 to be chosen too large. In such cases it may be 
necessary to restrict the maximum value that 6 can attain. That is, 
the method would be allowed to be self-adapting within a range of 

values between and some 6 . Further investigations of this 

± max 

scheme on various problems are recommended in order to evaluate its 
effectiveness as a general approach. 



Numerical Investigations With Distinctive 
Features of the PSPM 

\ 

Besides the convergence modifications, there are several distinc- 
tive features of the PSPM that may cause it to operate differently from 

t 

other Perturbation and Quasilinearization methods. One of these fea- 
tures is the capability for forcing the solution of the linearized 
equations to satisfy given terminal constraint functions to a specified 
accuracy at each iteration. This capability was used only in the ter- 
minal stages of convergence for the results presented and had no pro- 
nounced effect on whether convergence was actually obtained. It was 
found that the most optimum use of the capability was to restrict the 
PSPM to use only one Newton-Raphson iteration when the initial condi- 
tion step size norm was greater than 6, and to use no more than two to 
four Newton-Raphson iterations when this norm was less than 6. By 
restricting the PSPM to use only one Newton-Raphson iteration at all 
PSPM iterations, the method was operated in a fashion very similar to 
the Perturbation methods discussed by Lastman [26] and Lewallen [2U]. 
The primary difference in the normal PSPM operation and the restricted 
operation was that one to three total trajectory iterations were saved 
in the normal operation mode at the expense of one to five extra 
Newton-Raphson iterations. It is believed that the fewer trajectory 
iterations required resulted from better final time estimation obtained 
during the last several iterations. A definite savings in computer 
time was realized since the computer time required for a Newton-Raphson 
iteration is small compared to the time required for a total trajectory 
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iteration. This savings in computer time averaged about 20 percent for 
the cases compared. This result was not consistent for all cases and 
there was some dependence on the initial value of 6, since this 
affected the point in the terminal convergence phase where extra Newton- 
Raphson iterations were started. 

The generality of the PSPM operation makes it possible to use 
initial guess values for the a's other than the 1, 0, 0 values dis- 
cussed previously. An examination of the dtu/dt^ terms of the 
Jacobian matrix of equation (3.7) reveals that when the 1, 0, 0 values 
are chosen, only the reference solution influences these elements of 
the matrix on the first Newton-Raphson iteration. However, the influ- 
ence of the perturbed particular solutions can be obtained on the first 
iteration by assigning "weighting factors" to the various solutions 
with the initial choice of the a's. A typical choice investigated for 
example problem 1 was = 0.4, a^ = 0.3, = 0.3, so that the sum 

of the values totaled to 1 and more "weight" was given to the reference 
solution n p . During terminal stages of convergence, the initial 
guess was switched back to a^ = 1, = 0, a^ = 0. The results with 

this type of operation are inconclusive. In a comparison with the 
normal a selection procedure on a set of four different cases, this 
operation produced convergence in fewer iterations for two of the cases 
and required more iterations for the other two. This unique feature of. 
the PSPM makes it more general than other Perturbation and Quasilinear- 
ization methods, and it may be found to be more useful for other 
problems . 



Another distinctive feature of the PSPM investigated was the capa- 
bility for using different perturbation factors, 8^ and appearing 

in equation (4.22). Results were compared on various starting vectors 
using all 8 = 1.2 in one case and all 8^ = 0.5 in the other. All 
y^ were selected to be 0.1. If only one Newton-Raphson iteration at 
each solution iteration of the PSPM is made with the standard 1, 0, 0 ... 
guess on the a's, then theoretically the results with different per- 
turbation factors should be identical. However, significant differences 
were noted due to purely numerical causes. • These differences were sig- 
nificant enough to cause the method to require a different number of 
iterations for convergence when different perturbation factors were 
used. However, this difference was never more than one or two itera- 
tions. The results do point out the importance of minimizing numerical 
round-off errors in the matrix inversion computations. In this connec- 
tion, an important advantage results from using the particular solution 
method for solving the linear system, since the user can exercise con- 
trol over the numerical values which form the Jacobian matrix in equa- 
tion (3.7) by selection of appropriate perturbation factors. 

Results With Example Problem 2 

The second example problem, having a higher dimensionality and more 
complex terminal boundary conditions, would appear to be a more diffi- 
cult problem to solve than the first example problem considered. How- 
ever, once the first example problem is solved, the difficulty of 
guessing Lagrange multipliers for the second example problem is greatly 
reduced. The family of problems obtained by considering optimal 
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trajectories for different launch dates is most easily parameterized by 

the value 8 , the central angle of Mars at the launch time t . For 
o .■ o 

the "open problem" discussed previously, the Lagrange multipliers and 

final time for A_(t ) = -1 were found to be 
3 o 

X l( t Q ) = -0.494865 
A 2 (t Q ) = -1.07855 

*4 = 0 

t f = 3.319437 

and the corresponding .value of 8 q is easily determined to be 
0 Q = 0.7264 radian by using the time of flight, the angular velocity 
of Mars, and the final central single of the spacecraft in the open 
problem. To solve the second example problem for any value of 0 , a 
succession of problems having initial Mars central angles defined by 

0 =0 +A6 i = 1, 2, 3 ... 

°i i-1 ° 

where A0 q is some small increment, is solved in sequence using the 

converged values- of the previous problem as starting guesses for the 

next. The process is continued until a solution with the desired value 

for 0 q is obtained. The convergence characteristics of the PSFM were 

investigated on this example problem by studying allowable magnitudes 

for A0 . 
o 

Using the converged values for the open problem, a solution was 
first obtained for 0 q = 0.8. Repeated attempts to solve the problem 
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with 0 q = 0.9 using guess values from the 8 q = 0.8. solutions ended 

in failure. It was decided that A0 =0.1 was too large and A0 

o , o 

was reduced to 0.02. After solving several problems with this increment 
for A0 q , difficulties were again encountered. After reducing A0 q 
further to 0.002, the sequence of converged problems shown in Table III 

TABLE III 

CONVERGED MULTIPLIERS AND FINAL TIME FOR 
VARIOUS INITIAL MARS LEAD ANGLES 


Lead angle 0 , 
rad 

X 1 

X 2 ' 

X 4 

*f 

0.726k 

-0.4948 

-1.078 

0.000 

3.3194 

0.800 

-0.2321 

-1.994 

- 0.5177 ' 

3.3586 

0.820 

-0.0638 

-2.58 

-0.8379 

3.3776 

0.840 

0.2434 

-3.64 

-1.4106 

3.3985 

0.860 

0.9824 

-6.1718 

-2.7660 

3.4208 

0.880 

5.1917 

-20.5144 

-10.4111 

3.4443 

0.882 

7.0329 

-26.7821 

-13.7496 

3.4467 

0.884 

10.5011 

-38.5871 

-20.0369 

3.4491 

0.886 

19.4653 

-69.0958 

-36.2848 

3.4515 

0.888 

• 96.8223 

-332.3584 

-176.4839 

3.4540 


was obtained. The data in the table relate 0 q with the corresponding 
converged values of multipliers and final time. An examination of 
these data indicates that the PSPM was displaying good convergence 
characteristics on the boundary value problem but was getting nowhere 
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with finding' a solution to the optimization problem defined by 
0 q = 0.9* After plotting the data in Table III versus 8 , and noting 
the asymptotic character of the Lagrange multipliers as 0 q approached 
0.9, it was realized that all multipliers were seeking large values 
with respect to the normalized value of -1 for A^Ct^. This suggested 
that the unnormalized value of A^(t Q ) approaches zero as 0 q 
approaches 0.9. Since it was known that A^ would not be zero for this 
problem because of the constrained final central angle 0(t^.), the 
multipliers were normalized to A^ = -1. This eliminated the diffi- 
culty with convergence. 

With the problem normalized to A^ = -1, it was found that a value 
of A8 q =0.5 radian could be used to generate optimal trajectories 
for 0 q = 1.0, 1.5, 2.0, ... 6.5, 7.0 with an average of 11 iterations 
per problem. In this study, a maximum step size norm, 6, equal to 0.5 
was used without the self-adapting feature previously discussed. Opti- 
mal trajectories for 0 q < 0.T26U were also obtained. In this case it 
was necessary to normalize the multipliers to A^ = +1 in order to 
obtain the proper sign relationships between the multipliers. A plot 
of converged multipliers and final time as a function of 0 q is given 
in Figure 7. 

The good convergence characteristics of the PSPM were also demon- 
strated for this example problem by solving the problem for 0 q = 3 
using initial guess values for A^, A^, A^, and t f from the con- 

verged values of the problem with 0 q = 1 in 10 iterations. The. 
difference in the initial guess . traj ectory and the final converged 
















81 


trajectory is illustrated in Figure 8, where several of the optimal 
trajectories for different initial values of 0 Q are displayed. 

Interesting features evident in Figure 8 are the two distinctive 
types of optimal trajectories which result from launching before and 
after the most favorable launch date which corresponds to the open • 
problem (6 = 0.126b radian). This behavior has been previously dis- 

cussed by Kelley [53], who solved this problem with different numerical 
values for thrust and initial mass using a direct optimization method. 
The trajectories corresponding to early (0 >' 0.7261+) launch dates have 

a "pursuit from behind" character, while the spacecraft when departing 

from late (0 < 0.7261+) launch dates tends to "wait" for Mars to over- 

o 

take it. The severe time-of-f light penalty associated with not launch- 
ing on the most favorable date is shown in Figure 7. It is also evident 
from Figure 7 that the "pursuit from behind" type of trajectory has a 
shorter transfer time than the "waiting" type for most of the unfavor- 


able launch dates. 





CHAPTER VI 


CONCLUSIONS AND RECOMMENDATIONS 

Several important extensions and modifications to existing indirect 
optimization methods have been made. The method of particular solutions 
was extended in order to solve linear boundary value problems with 
boundary conditions specified in the form of nonlinear functions of the 
dependent and independent variables . This extension was incorporated 
into a Perturbation method for solving the nonlinear boundary value 
problem which results from formulating optimal control problems for 
solution by an indirect method. This new Perturbation method, called 
the Particular Solution Perturbation Method (PSPM), reveals a new 
approach for solving problems with unknown final time which can reduce 

I 

the number of trajectory iterations required for convergence to the 
optimal trajectory. The application of this new method for treating 
unspecified final time problems was simplified by the use of a power 
series numerical integration method which was ideally suited for the 
forward and backward variable step integration required. The method is 

not restricted for use with power series integration, however, and may 

/ 

be implemented with any numerical integration scheme. 

The PSPM was found to have excellent convergence characteristics. 
The range of convergence of the indirect optimization approach was 
extended far beyond that of previous methods without compromising the 
rapid convergence of this approach, and thus now places the indirect 
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optimization approach in a more competitive position with direct 
methods. Although the PSPM utilizes several features not available 
with currently more well known indirect optimization methods, the ex- 
cellent convergence characteristics are primarily due to the easily 
applied modifications of equations (4.13) to (4.l6) together with upper 
and lower bounds placed on allowable values of the unknown final time. 
Thus, it is expected that with these modifications, other indirect 
optimization methods currently programed need not employ the particular 
solution method presented and the more unfamiliar power series integra- 
tion in order to ob'tain the good convergence characteristics displayed 
by the PSPM. 

As a result of this study, several areas are recommended for 
future investigations. Although the good convergence characteristics 
of the PSPM are not believed to be unique to the example problems pre- 
sented, the convergence characteristics of this method should be studied 
on other problems of larger dimension and of a different nature (such as 
atmospheric reentry problems with inequality constraints on control and 
state variables) in order to support the claims made here. 

It would appear that the use of regularizing transformations 
discussed by Tapley, Szebehely, and Lewallen [^4] would be as beneficial 
with power series integration as with more conventional integration 
schemes in solving trajectory optimization problems. This should be 
investigated. 

The methods presented here for solving two-point boundary value 
problems are not restricted to the typical problem which results from 
trajectory optimization. The modifications employed to extend the 
range of convergence should be equally as beneficial on any two-point 
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or multi-point boundary value problem solved by a Perturbation or 
Quasilinearization method, and this should be investigated. The advan- 
tages of power series integration may be more fully realized for multi- 
point boundary, value problems because of the ease with which the method 
can obtain solution values at any value of the independent variable. 

The details of adapting the method for solving nonlinear multipoint 
boundary value problems have been worked out in this study, and several 
such problems should be solved to test the usefulness of the power series 
method. The Newt on-Raph son method utilized for solving unspecified 
final time problems could be applied to multipoint boundary value 
problems where several boundary conditions at unspecified values of the 
independent variable are known. This should be demonstrated. 

Finally, this investigation has revealed the Perturbation approach 
to have several practical advantages over the Quasilinearization approach 
for solving nonlinear boundary value problems. Besides requiring fewer 
integrations per iteration and less computer storage than the Quasi- 
linearization method, the Perturbation approach admits the capability 
for simultaneous integration of the reference solution and linearized 
equations, which in turn allows for variable step integration and 
capability for extreme solution accuracy. However, the convergence of 
the Quasilinearization approach has been rigorously established [ 19 ], 

[ 29 ] , [ 34 ] for boundary value problems of a less general nature than 
considered in this investigation, while the Perturbation approach is 
lacking in this regard. When compared in numerical studies [ 24 ], [ 42 ], 
the methods have displayed, similar convergence characteristics, and the 
Perturbation approach as modified in this study exhibits convergence 
characteristics superior to the Quasilinearization method reported in 



reference [35] for the same example problem. Further theoretical in- 
vestigations of the Perturbation method are needed to establish the 
necessary and sufficiency theorems for convergence which must exist. 

In the past, theoretical investigations of the Quasilinearization method 
have been easier because the method involves iterative solutions of a 
system of linear differential equations only, while the Perturbation 
approach involves iterative solution of both linear and nonlinear 
systems. However, it was established in this investigation that the 
nonlinear solution at each iteration of the Perturbation method is a 
particular solution of the linear system. In addition, as shown in 
Appendix C, initial and final values of the nonlinear reference solution 
can be related through the same fundamental set of solutions used to 
construct the general solution of the linear system. Perhaps some 
advantage can be made of these properties in future theoretical in- 
vestigations of the Perturbation method.’. 
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APPENDIX A 


REDUCTION OF. AN OPTIMIZATION PROBLEM TO A TWO-POINT 
BOUNDARY VALUE PROBLEM 


For the example problem considered in Chapter V, it is necessary 
to determine the optimal thrust vector control for a constant low- 
thrust rocket in a planar Earth-Mars orbit transfer so that the trans- 
fer is completed in minimum time. The orbits of both Earth and Mars 
are assumed to be circular in this example. In this appendix, the 
necessary conditions for optimal control outlined in Chapter II are 
applied to the example problem considered in Chapter V in order to 
reduce the optimization problem to a two-point boundary value problem. 

The equations of motion for the thrusting rocket , expressed in 
a polar coordinate system with origin at the sun, are given by: 


u 





r m 


r 


u 



m 


-c 


where T is the thrust magnitude, GM is the solar gravitational constant 
c is the constant mass flow rate of the rocket exhaust, and B is the 
time varying thrust control angle. 
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In order to apply the necessary conditions outlined in Chapter II, 
the following substitutions are made: 



So that the equations of motion are written in the form x = f(x,u,t) 
corresponding to equation (2.1), 

2 

X 2 GM T . . , 

— - — r + — sin u.. = f. U,u,t) 

x_ 2 x c 1 1 

3 x 3 5 

• — X 1 X 2 T 

x 0 = + — - cos u.. = f 0 (x,u,t) 

2 x^ x^ 1 2 

x 3 = x ± = f (x,u,t) (A.l) 

. x 2 

X IT = ^* f l. (x ' u - t) 

x 5 = -c = f (x,u,t) 

The Hamiltonian function is given by 




9 ^ 


and the Euler-Lagrange equations are obtained from the necessary 
conditions, equation (2.5), 

/ \T 
• _ (3H\ 

which, after some simplification, can be written 


X. = 




X„ = 


~\3 / 1 *XsJ 2 " i X 3/ 


X. = 


X, = 0 


2GM 

3 


x l x 2\ / 2\ 

X l-(— X 2 + “K 

X 3 / \3/ 


x 5 = - — ^X ± sin u l + X 2 cos u^ 


The control variable u^ is eliminated from equations (A.l) 
(A. 2) by application of necessary conditions (2.6), 


3H 

3u n 


= 0 


with the result 


Xj — cos u_, \ 

\ X 5 


- X 2 fe sin 




(A. 2) 


and 


0 



Simplifying , 
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- — = tan u 
A 2 1 


which implies , 


sm u, = 


1 


cos u., = 


1 4 


2 2 
X 1 + X 2 


(A. 3) 


The necessary Weierstrass condition 


E = F(t,x,X,U,X) - F(t,x,x,u,X) - — (t ,x ,x ,u,X ) (X - x) 

i 3x 


- — (t,x,x,u,X)(U - u) > 0 


J 1 1 | ^ 

where F = X (f(xjU,t) - x) and X and U ( are nonoptimal but per- 
missible values for x and u, is imposed to resolve the ambiguity in 
sign appearing in equations (A. 3). Since the equations of motion 
must be satisfied on a permissible trajectory, 

F(t,x,X,U,X) - F(t,x,x,u,X) = 0 

and since, from the optimality condition, 


3F _ 3H 
3u 3u U 



the Weierstrass E-function reduces to 


96 


E = - (t ,x ,x ,u,A ) (X - x) > 0 
3x 

which for this problem simplifies to 

T, • • v 

E = A (X - x) > 0 

Substituting with x = f(x,u,t) and X = f(x,U,t) 


[t , . ~ 


"t . Tt . 

— /sin U_ - sm u_\ 

+ A 

— /cos U, - cos u.,\ 

H 

— 1 

ir 

X 

2 

X { 1 1) 



y 


Substituting for sin u^ . and cos u^ with equations (A. 3) yields, after 
some manipulation. 





If the above expression is to be nonnegative for all admissible 
values of U^, then the negative sign on the radical must be chosen, 
hence , 


sin u^ 


cos u^ 



Substituting the above expressions into equations (A.l) and (A. 2) 
eliminates the control parameter u^ from the equations of motion. 
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The terminal boundary conditions , ij>[x(tp ,t^] , which correspond 
to necessary conditions (2.8), are given by 


. h * - 0 ■ 0 

*2 = X 2( t f) - V M = 0 
* 3 = x 3 (t f ) - r M = 0 

V M 

*k * *u(*t) - 9 M (‘f) - x u(M - 9 m(‘o) - J- ‘f ■ 0 


where the subscript M refers to the value for Mars. The performance 
index <t>[x(t^,) ,t^] ' "is simply t^ for a minimum time transfer. 
Applying necessary condition (2.7), 


yields 



u i - hfr) = 0 

U 2 - X 2(‘f) * 0 
w 3 - X 3 (*f) ’ 0 

- x l .(h) * 0 

- S(‘f) - 0 
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One version of the example problem in Chapter V assumes that 
the final angle x^t^ is not constrained. In this case it is 
seen from the above application that 

H^f) vH ' 0 

just as A ^t ^ was determined to be zero since no constraint was 
placed on the final spacecraft mass. 

The necessary condition corresponding to equation (2.9) 


becomes 


1 + v 


+ Vl^f) * x 5 i 5'( t *) = 0 

Using equations (A. 5 ) 5 can be eliminated from the above equation 

to yield 


, M . 

1 + ♦ » 2 X 2 + X 3 x 3 * - — j + X 5 x 5 


= 0 
t , 


(A. 6) 


-* f 


Since the constant Lagrange multipliers v have been eliminated from 
the terminal boundary conditions, there -is no reason to compute them 
and the trivial differential equations v = 0 appearing in equa- 
tion (2.10) can be eliminated from the formulation of the two-point 


boundary value problem. 
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Equations (A.l) and (A. 2) provide 10 ordinary differential equations 
to be solved with five initial conditions provided by the known initial 

position, velocity, and mass of the spacecraft as it leaves an assumed 

/ 

circular Earth orbit about the sun. 

*1 (‘o) = 0 
x 2 (*" o] V Earth 
x 3 (^o) r Earth 

' ■ . x -(\) ■' 0 ' • 

x 5 ( t o) = "o 

Equations (A. 5) yield a zero value at the final time for X^(t f ), 

Since t f is not specified, five additional boundary values are re- 
quired for the 10 differential equations. Four of these conditions are 
provided by equations (A. 4) and the fifth condition is provided by 
equation (A. 6). 

From a computational point of view, it is desirable to normalize 
the values of the variables in the differential equations so that some 
degree of numerical magnitude compatibility is achieved. Since it was 

desired to compare. the. numerical, results of this .investigation. with 

previously published results of reference [24], the normalization scheme 
of reference [24] was employed. In this scheme, the fifth equation in 
(A. 2) is eliminated together with terminal condition (A. 6). The 
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initial values of the Lagrange multipliers are normalized to the 

unknown value of A_(t ) such that X _ ( t ) is specified to be -1. 

3 o 3 o 

This is possible since the remaining Euler-Lagrange equations (A. 2) 
are linear and homogeneous in the A’s, and since only the ratios of 
A^ and A^ appear in the equations of motion (A.l). In addition to 
providing better numerical accuracy, this normalization technique also 
reduces the complexity of the boundary value problem since an additional 
initial condition is obtained and the terminal condition (A. 6) need 
not be used as a boundary condition. Since equation (A. 6) must be 
satisfied on the optimal trajectory, it can be used to recover the un- 
normalized values of the Lagrange multipliers. However, there appears to 
be no practical reason to recover these values. Other normalized values 
of parameters of interest are 

Gravitational constant of the Sun, GM = 1.0 
Initial spacecraft mass, m o = 1.0 

Initial spacecraft velocity, v g ar1 -^ = 1-0 

Initial spacecraft radius, r„ ,, =1.0 

Earth 

Terminal spacecraft velocity, v ^ ars = 0.81012728 

Terminal spacecraft radius, r ^ ars = 1-5236790 

\ 

Thrust = 0.14012969 

Mass flow rate = 0.074800391 



With these normalized values, units of various physical quantities 

are 

. Length unit = 1 astronomical unit = 1.495987 x 10 11 meters 

Mass unit = 6^7978852 x 10^ kilograms 

4 

Velocity unit = 2.9784901 x 10 m/sec 
Force unit = 4.0312370 newtons 

Time unit = 5-0226355 * 10 6 second = 58.132355 days 



APPENDIX B 


A POWER SERIES NUMERICAL INTEGRATION METHOD 

One of the most important facets of any method for solving multipoint 
boundary value problems is the numerical integration scheme used. Since 
the numerical integration of differential equations consumes the bulk of 
computer time, it is desirable to have a fast and accurate integration 
method. Among the most popular integration methods are the well-known 
Runge-Kutta formulas and predictor-corrector methods. In this appendix, 
a power series integration method is presented which has several features 
which make it uniquely suited for use as an integration method in solving 
multipoint boundary value problems. - 

The capability for solving differential equations by power series 
expansions has been known since B. - Taylor (1685-1731). However, this 
method has not enjoyed the popularity of other numerical integration 
schemes. This is probably due to the fact that it is an impractical 
method 'for hand calculation or even desk calculators, and thus did not 
receive the early attention and development of the currently more 
popular integration methods used on digital computers. 

With modern digital computers, the cumbersome application of the 
power series method is easily overcome and its practicality is evidenced 
by its high accuracy, large step sizes, good speed, and variable step 
size capability. The use of power series as a method for digital 
computers has been studied by Fehlberg [55] and Hartwell [ 5 6]. Detailed 
programing steps for the method are given by Doiron [57]. The method is 
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reported to have superior speed and accuracy characteristics on problems 
where integration must be made over large intervals of the independent 
variable with frequent changes in the integration step size. 

Given a system of differential equations , 

z i " 

1=1. 2. ••••< (B.D 

the method assumes a power series expansion exists in a neighborhood of 

t for each of the variables z. of the form 
o '1 


? i (t) Z 1 <1C> ' t o) k ' 1 (B - 2> 

k=l 

(k) 

where the z^ are power series coefficients and t Q is the value 
of the independent variable where the power series expansion is made. 

In the following, it will be assumed that power series solutions of 
equation (B.l) exist. 

Letting (z^) denote the power series coefficients of z^, 
it is easily determined from term-by-term differentiation of (B. 2) that 



which yields a recurrence relation for (z^ )^ +1 ^ if (z^)^ is 
known. Since z^ = *x^ Z ±' Z 2' ‘ * Z N’^ * ^ ^°^ ows that power series 
expansions of the functions f^ exist with power series coefficients 



denoted by f . Thus, by equality of power series, we have equality 
of the power series coefficients, and (B.3) can be written 


( 2 0 


(k+1) _ _£i) 


(k) 


(B.4) 


The main difficulty in applying the method is involved with deter- 

(k) 

mining the coefficients f . The functional relationship between 

the -z , as expressed by (B.l), must be carried out in "power series 

arithmetic" and when coefficients of like powers of (t - t ) are 

o 

(k) 

collected, these coefficients represent the f^ . It can be shown 

(k) 

(see Theorem 13-27 of Apostol [ 58 ]). that the coefficients f^ involve 

(k ) 

only the first k coefficients, , of the power series for the z^. 

This guarantees that equation (B.4) does actually represent a recurrence 

(k+l) (k) 

formula for z^ in terms of the first k, z^ . The application 

of the power series arithmetic is greatly simplified through the use of 
auxiliary series and repetitive application of known algorithms for 
series addition, subtraction, multiplication, and division. Easily 
programed algorithms are also known for generating power series coeffi- 
cients of transcendental functions of power series such as sin z and 
B 

z , where . 3 .is some real number. 
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Let u, v, and w be power series of the form 


u -£\( t - t o) 


k=l 


=£\(‘- t o) 


k-l 


k-1 


k=l 


“ “"( t ' *«) 


k-l 


k=l 


The following power series operations are defined by: 


Addition: 


w = u + v=^w = 


Subtraction: w = u - v =^w = 


Multiplication: w = u • v =^w = 


\ + v k 

\ " v k 
k 


u.v, 

1 k-i+1 


Division: 


w = u/v 


L-l 


W k = 


\ - 7 =f w i\-i + i 


U-, 

» W = 1 


Square root: 


” = <^=* U k = 55^ ^ Vk-i+lj 


"l ’^1 


r 



io6 


Integral or fractional powers: 

6 » 6 
W = U =? W^ = 


w 2 = 8 u 2 w 1 


for k > 1 


w. 


k+1 ku. 


k-1 


st Vi'i * £ li(B + 1 


i=l 


> - ‘Vk-W 


Sine and cosine: 


u = cos w u, = cos w. 


v = 


sin w v. 


= sin w 


Vi 


1 

1 k 

k E iw i + iv i+ i 

i=l 


k+1 


1 k 

k Xw iW i+l U k-i+l 
i=l 


These algorithms are sufficient for the differential equations which are 


solved in this thesis. 
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An Example Problem . 

Given the differential equations. 


x = y + 


cos w 


(* 2 ♦ 7 ) 


3/2 


y = x + 


sin w 


C 2 ♦ /) 


3/2 


w = w * y 


with initial conditions, x(t Q ) = x ^» y(t Q ) = y^, and w(t Q ) = w^, 

! 

generate power series coefficients for x, y, arid w for a power 

series expansion about t . ! 

o 

The coefficients are obtained by computing in sequence the k 

coefficient of each of the auxiliary series in equations ( B . 5 ) using the 

\ . 

above algorithms, and then applying the recurrence relations (B.6) to 
obtain the (k + l)st coefficients for x, y, and w. The process 
is repeated for k = 1, 2, 3 ... N-l, where N is the desired number 
of coefficients. j 


a = 

w • 

y 

\ 

u = 

cos 

w 


v = 

sin 

w 


b = 

x • 

X 


- - 



- - -1 


c = y • y ) 


r = b + c 



e = u • s 
f = v • s 


(B.5) 
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Vl = 


M - 4 m 





(B.6) 


Once the desired number of coefficients N are computed, the next 
step in the integration process is to determine the integration step 
size (t - t Q ), which can be used with the available coefficients, while 
maintaining a specified numerical accuracy in the evaluation of the power 
series solutions. In general, a larger step size may be used with a 
larger number of coefficients. A practical limit for the number of 
coefficients which should be computed is determined by the magnitude of 
the N coefficient of the series to be evaluated. Digital computers 
have a largest and smallest value of the magnitude of a number which can 
be accurately represented. Depending on the radius of convergence of 
the series, the coefficients may approach one or the other of these 
limits. The number of coefficients computed should be limited to avoid 
these number magnitude limits. 

It is assumed that N power series coefficients are available for 
evaluation. A method for determining the largest allowable integration 
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step size is developed below. The equivalence of the power series 
expansion and Taylor series expansion of a function x(tj about a 
point t is well known. The Taylor series expansion can be written 


xU) = x(t o ) + x'(t o )(t - t o ) + %(t 0 ) (t - t o ) 2 * 

+ ( l - l o) N * >W 

where the remainder, is given by 


H+J. /. v 

x ( l) /. . X N+l 

.. ’Sl+l * (n+l)l (* - * 0 ) * *1 •= 11 

If it is assumed that x N+ ^(t^) ^ x^ + ^"(t ), then the truncation 

error is of the order of the first neglected term in the 

summation of the series. Since in most applications it is desirable to 

limit the relative truncation error rather than the absolute truncation 

error, the step size (t - t ) should be chosen to meet a specified 

relative truncation error bound, e , . Let the relative truncation 

rel 

error be . defined by 


'rel 


^numerical X ^exact 


x(t) 


exact 


For all practical purposes, the denominator in the above expression 


need only represent the magnitude of x(t), and, therefore, it may be 



no 


approximated by the summation of the first p terms of the power series 
representation of x(t). In light of the foregoing assumptions, if 
x(t ) is to be accurately represented by 

x(t) \( t - ‘of' 1 

k=l 


then a reasonable truncation error check would require that each of the 
last several terms (r “terms, for example) of the summation be less in 
magnitude than 

s* E - ‘of 1 

k=l 

where e rel is the specified relative error allowable. The value of 
r depends on the severity of the test desired. In practice, r = 2 
has been sufficient to maintain the desired accuracy. With r specified, 
a requirement of the test can be stated. 


*N-r+l 


(‘ - ‘o) 


N-r 


< e 


rel 


N-r 

£ ’ ‘°) 
k=l 


k-1 


(B.7) 


It is desirable to solve for (t - t Q ) which will satisfy this 
test. Since the summation on the right of the inequality (B.7) is used 
only to approximate the magnitude of x(t), let the magnitude be 
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approximated by the constant term in the expansion x^. Then, taking 
logarithms of both sides above yields 



Exponentiating both sides yields an estimate for the allowable convergence 
interval 



lo s( e rel 

X 1 

\" 

*N-r+l 

) 

(N 

-r) 



(B.8) 


Since normally there are more than one series to be evaluated!, it 
is necessary to determine the largest convergence interval common to all 
series to be evaluated. Noting in (B.8) that |t - t Q | is a monotonic 
increasing function of | x 1 / x u_ r+1 1 » it is only necessary to compute 
for each series the value analogous to |x^/Xjj r+ ^| determine which 
of the series will yield the smallest convergence interval. A trial 
convergence interval can then be obtained by multiplying the quantity 
on the right of (B.8) by a positive number, P, less than unity to insure 
the inequality. A trial step size determined in this manner may still 
not satisfy the convergence test ( B . 7 ) because the approximation of 
x 1 for 



N-r 

EM* 


k=l 



may not be sufficiently accurate. Since a failure of the test (B.7) 
during evaluation for any one series would require reducing the convergence 
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interval and reevaluating all series, it is desirable to choose the 
number P so that the likelihood of convergence test failure is very 
small. The choice of P is dependent on the number of coefficients 
used, since longer convergence intervals allow for a larger change in 
magnitude of the solutions. In practice, values of P of 0.9 or 
0.8 have worked well with 10 to 20 term power series. It should be 
noted that special logic is required in the case where either of the 
coefficients x^ or x^ . is zero. 

Returning to the example problem, once a trial convergence interval 
At for the series x(t), y(t), and w(t) has been determined, the 

series can be evaluated for any value t^, |t^ - t Q | < At. If the 
solution is required for some t^ outside the convergence interval, the 
series may be evaluated at t^ = t Q ±At and new series expansions about 
tg can be obtained. The analytical continuation can be repeated until 
power series which will converge at the desired value of t are obtained. 



APPENDIX C 


EQUIVALENCE OF TWO METHODS FOR MODIFYING BOUNDARY 
CONDITIONS OF LINEAR DIFFERENTIAL EQUATIONS 

The following derivation is made to support the claim made in 
Chapter IV regarding the equivalence of two methods for determining 
modified initial conditions for the reference solution of the Particular 
Solution Perturbation Method. To simplify notation in the derivation, 
the general solution of the N dimensional linear system (U.12) with 
m specified initial conditions is written as a linear combination of 
particular solutions 

S+i 

n y(t) = Ect _P (t) = P(t)a S = N - m 
k=l * n n 

where n P(t) is an N by (S+l) matrix formed with columns of 
particular solutions ^ 

P(t) = [ p 1 p 2 ... p S+1 l 

n [rr n n J 


and a is the vector of superposition constants 
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Ilk 


Let terminal boundary conditions at the fixed final time be specified 


Hy i ^f) =Z fi i 8 k + 1, k + 2, ...k + S 

for some 0 < k < N - 1. An (S+l) vector y(t) and (S+l) by 

(S+l) submatrix [p(t)] are formed by partitioning out the (k+l)st 

s 

through the (k+S) rows of n y(t) and n P(t) and then augmenting each 
by a row of l's to obtain 


y(t) 


1 


1 1 , ... 1 . 

n 


1 . 2 S+l 

y k+l 


n P k+l n P k+l n P k+l 

n 

[P(t)] Q = 

1 2 

y k+2 

• 

• 

• 

s 

n P k+2 n P k+2 

• 

• • 

• • 

• • 

n 

y k+S 


S+l 

n F k+S n y k+S n p k+S 


Method 1 


Given terminal boundary conditions z ^». i = k + 1, ... (k + N - m) 
written as an S = N - m vector z^, modified terminal boundary condi- 
tions for the system (4.12) are formed by 



ez f + (1 - e) n z f 


(C.l) 
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where n y f and n z f are S vectors with elements 


y k+l^ t f ^ 


\n<V 

"WV 


Z k+2^P 

• 

• 

and 

• 

• 

• 

°W t f ) 


• 

"WV 



- - 


respectively, and n z(t) is the nth reference solution of 
method. 

It follows from applying the boundary condition (C.l) 
S+l 

auxiliary conditions o, = 1 that 

k=l 


Wt) = [ p (*f)] “ 

S 



Solving for a , 



a Perturbation 


and the 



ll 6 

Initial conditions for the (n+l)st reference solution are obtained from 


n+1 


4o) = Ho) ■ n P ( t o) a - M M 


-1 


n 


'fj 


(C.2) 


Before proceeding with an analysis of the alternate method for 
computing z(t Q ), it is necessary to establish an identity which will 
be needed. Using Theorem 1 of Chapter IV, it is possible to write 

n z(t) = ^(t) 

subject to 

Ho) “ /( l o) 

Forming an (S+l) vector z(t) by partitioning n z(t) and then 

augmenting with a 1 in the same manner that y(t) was formed from 

n y(t) yields z(t) « [P(t ) ] y, where y = [1,0,0,. . .0]^. At the final 

s 

time 

= PH]; ■ 
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so that 



which allows initial 


and final conditions for 




to he related by 




n P &> * [n^o) [ P (‘f3 



(C.3) 


This seemingly awkward result will allow considerable simplification in 

/ 

the derivation which follows. 


Method 2 

For this method, n y(t Q ) is computed without modifying terminal 
boundary conditions, and then z(t Q ) is confuted by 



n 


z 


(tj ♦ t( n y(t o ) 



(C.U) 


It is necessary to show that z(t ) computed in this manner is equal 
to the result, (C.2). 

First, n y(t Q ) is determined by applying the unmodified terminal 
boundary conditions. Using similar notation as before 

rii 


y (V ■ ( P (M 


s 
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so that 


-1 

r 1 

[ p ( tp )] s 

_ z f_ 


and consequently 


"ifo) ■ n P ( t o)“ - [n P ( t o)j [ P ( t f2 ' 


I f. 


Implementing the modification (C.H) 


n+i 


= “*(0 + e HvL 


*- y- 


- MOi 


= (1 - £)"*( t c ) * =[ n P(t c )] 


-l 

s 


Using the result (C.3) and factoring 


n+1 


*W = [n P (‘o1 [ P ( t f)]. 1 * 


'fj 


+ (l - e) 


n 

L z fJ 


Substituting with equation (C.l), the final result is obtained 


n+1 


W ■ [/W] [ p ( t fl * 


(C.5) 
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A comparison of equations (C.2) and (C.5) reveals that they are identical 
and therefore Methods 1 and 2 are shown to be equivalent. 
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